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Abstract 

The activities presented are a broad based approach to advancing key hydrogen related 
technologies in areas such as fuel cells, hydrogen production, and distributed sensors for 
hydrogen-leak detection, laser instrumentation for hydrogen-leak detection, and cryogenic 
transport and storage. Presented are the results from research projects, education and outreach 
activities, system, and trade studies. The work will aid in advancing the state-of-the-art for 
several critical technologies related to the implementation of a hydrogen infrastructure. Activities 
conducted are relevant to a number of propulsion and power systems for terrestrial, aeronautics 
and aerospace applications. 

Hydrogen storage and in-space hydrogen transport research focused on developing and 
verifying design concepts for efficient, safe, lightweight liquid hydrogen cryogenic storage 
systems. Research into hydrogen production had a specific goal of further advancing proton 
conducting membrane technology in the laboratory at a larger scale. System and process trade 
studies evaluated the proton conducting membrane technology, specifically, scale-up issues. 
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1. Hydrogen Production Using Advanced Protonic Conductor 

Task PI: Dr. Eric Wachsman, Material Science & Engineering 

Gradute Students: Dr. Heesung Yoon (Post doc), Ta keun Oh, and Jianlin Li, Material Science 
& Engineering 

Research Period: August 3, 2004 to March 31, 2008 

Abstract 

A ceramic hydrogen (H 2 ) permeation cell was developed based on the advanced protonic 
conductor: europium-doped strontium cerate (SrCei-xEu x 0 3 . 5 ). The permeation cell is has a 
tubular design (with one end capped) to alleviate gas leakage issues common in planar designs. 
To maximize hydrogen permeation flux the basic cell architecture consists of a dense thin film 
SrCe^EuxOs-s membrane (for the actual hydrogen permeation) supported by a porous (for 
unimpeded gas flow) Ni-SrCe0 3 composite to provide mechanical strength to the thin film 
membrane and catalytic activity to the feed gases (from the Ni component). With this membrane 
cell a maximum hydrogen flux of 6.9 cm 3 /min was achieved. 

To optimize the performance of the membrane, the transport phenomena of SrCe^EUxOs-s was 
investigated, in dry and wet hydrogen atmospheres, as a function of temperature. To wit, the 
optimal europium dopant concentration in SrCe^EUxCh-s for maximum hydrogen permeation 
was determined from measurements of ambipolar conductivity, total conductivity and open 
circuit potential. Of the compositions of SrCe^xEuxO^g (x = 0.1, 0.15, and 0.2) studied, 
SrCe 0 .9Euo.i0 3 .5 showed highest total conductivity between 500 °C and 900 °C under both dry 
and wet hydrogen conditions. However, the highest ambipolar conductivity was obtained over 
the compositional range from SrCe 0 . 85 Eu 0 .i 5 O 3 . s to SrCe 0 . 8 Eu 0 . 2 O 3 . s , depending on temperature. 
Therefore, since H 2 flux is proportional to ambipolar conductivity these compositions were 
characterized for H 2 permeation efficiency. 

Finally, the use of the permeation cell to enhance hydrogen production from both water-gas shift 
(WGS) reaction and steam reforming (SR) was explored. Our SrCe 1 . x Eu x 0 3 . s -based permeation 
cell was employed to continuously remove the H 2 produced by the WGS or SR reactions from 
either reaction stream. The CO conversion (for WGS) and the methane conversion (for SR) as 
well as the total hydrogen production (for both WGS and SR) were considerably improved. 
Moreover, the use of the H 2 permeation cell allows the WGS reaction to be operated with low 
H 2 0/CO ratios, without the thermodynamic constraint. 

This technology will dramatically improve H 2 production technology from a broad array of 
conventional (natural gas, coal) and renewable (biomass) fuels. Moreover, it is a compact 
efficient technology that can be used to produce pure H 2 for aerospace transportation 
applications (jet fuel/JP8) once sulfur issues are dealt with. Finally, as this is a transition project 
it directly lends itself to technology transfer and commercialization. 
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Introduction 


Background 

Hydrogen is the most abundant chemical-energy resource in the world. Historically, hydrogen 
has been produced from fossil fuels, and for the foreseeable future, will continue to be so. The 
major source of H 2 is steam reformation of natural gas. Therefore, improvements in the 
efficiency and cost of H 2 production from natural gas are necessary in the near future. Gas 
separation membranes and membrane reactors, based on ion conducting ceramics, may 
provide the technological advance necessary to increase the efficiency and reduce the cost of 
H 2 production from natural gas. However, other sources of H 2 must be developed for the 
envisioned hydrogen economy, and coal provides the greatest U.S. domestic resource-based 
option. The U.S. DOE is developing a FutureGen plant based on coal gasification, solid oxide 
fuel cells (SOFCs), biomass resources and ion conducting membranes that will produce H 2 and 
electricity with zero emissions and carbon sequestration; thereby, not contributing to global 
warming. 

Membrane reactor technology holds the promise to circumvent thermodynamic equilibrium 
limitations by in situ removal of product species, resulting in improved chemical yields. Mixed- 
conducting oxide-membrane technology presents the possibility for a dramatic reduction in the 
cost of converting petroleum and coal derived feed stocks to hydrogen and other value-added 
hydrocarbons. The membranes are based on metal oxides (e.g., SrCe-i_ x M x 0 3 _ 8 , where M is a 
metal dopant, e.g., Eu) that exhibit both ionic and electronic (mixed) conductivity. Because of 
their significant electronic conductivity, these mixed ionic-electronic conductors (MIECs) have 
an internal electrical short and the ionic species selectively permeates through a dense film of 
the material under a differential partial pressure, as shown conceptually for H 2 permeation. The 
potential permeation flux rates of these materials are extremely high. 

High-temperature proton conducting ceramics are promising materials for many applications, 
including solid-oxide fuel cells, solid-state gas sensors, and membranes for hydrogen gas 
separation, etc. [1-3]. Recently, a series of high-temperature proton conducting ceramics, 
typically doped with multivalent cations, have been studied, based on BaCe0 3 . s [4-6], SrCe0 3 . 
s [7-9], and complex perovskites in the form of A 2 B'B''0 6 and A 3 B'B" 2 0 9 [10-12], BaCe0 3 . s based 
oxides exhibit oxygen ion conductivity comparable to their proton conductivity and thus are not 
desirable for hydrogen separation applications where oxygen is present [13, 14]. According to 
previous reports by Knight and Bonanos et a/., distorted orthorhombic structure of SrCe0 3 .§ 
inhibits oxygen ion conduction, thereby enabling high transference number for proton 
conduction [15, 16]. Therefore, SrCe0 3 based oxides have promise for selective hydrogen 
separation if their electronic conductivity can be improved by proper doping. 

The advancement of hydrogen separation technology has become important in energy 
applications for production of petrochemicals and hydrogen for fuel cells [2]. High temperature 
proton conductors based on perovskite-structured oxides have received increasing attention for 
hydrogen separation technologies. The separation process of protonic conductors is purely ionic 
transport, and results in nearly 100% pure hydrogen. Recently, our group studied Eu-doped 
SrCe0 3 for use as hydrogen separation membranes due to its mixed protonic-electronic 
conductivity [7, 8, 17-21], Eu-doped SrCe0 3 is especially promising since it provides high 
selectivity to protons. Also, T. Norby and Y. Larring showed the flux of each carrier is inversely 
proportional to membrane thickness [22]. Therefore, our research has been directed towards 
the development of thin film hydrogen membranes using tubular-type porous supports for 
increasing the hydrogen production by reducing the thickness of the membrane. Figure 1 
illustrates the basic design for producing pure hydrogen through the methane steam reforming 
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process using thin film proton conducting membranes. For a hydrogen production cell, a 
hydrogen permeable thin film is coated on the inner-side of a catalytic tubular-type porous 
support. A methane and steam mixture is flowed on the outer side of a reactor cell. At operating 
temperatures, CO, C0 2 , and H 2 are formed, and due to the hydrogen partial pressure gradient, 
H 2 permeates through the hydrogen permeable thin film. 




Membrane reactors can be used to produce hydrogen from hydrocarbon feed-stocks, from 
natural gas to coal to a variety of hydrocarbon containing waste streams such as landfill gas and 
swamp gas. The overall objective of this project is to demonstrate the feasibility of producing 
hydrogen from hydrocarbon based fuels using advanced proton conducting membranes. The 
objectives of this phase are to develop thin film proton conducting membranes on porous 
supports and to advance our fundamental understanding of these materials. The anticipated 
result is demonstration of high hydrogen fluxes through these thin supported membranes. 

Technical Approach 

The objectives of this report are to optimize electrical properties, to improve the stability of 
SrCe0 3 , and to fabricate and demonstrate a membrane reactor that produces hydrogen from 
hydrocarbon based fuels (natural gas, jet fuel, coal gas, or biomass) using advanced proton 
conducting membrane materials. 

Optimization of electrical properties for hydrogen permeation-Jo find optimized dopant 
concentration in SrCe0 3 , we conducted a study of the electrical properties of SrCe-i-xEu x 03. 8 (x = 
0.1, 0.15, and 0.2) in various oxygen and hydrogen atmosphere. The protonic conductivity was 
determined by impedance spectroscopy which allows one to extract bulk resistance from 
complex impedance as a function of frequency. 

Improvement of chemical stability ofSrCe0 3 - Zirconates, such as SrZr0 3 , show good 
chemical stability and mechanical strength and they are more stable against carbon dioxide gas, 
which reacts with cerate materials below 800 °C [23, 24], though the conductivity of zirconates 
is lower than that of the cerates [25], However, it is possible to achieve a compromise between 
the proton conductivity and chemical stability using a mixed solid solution of SrCe0 3 and SrZr0 3 . 

Fabrication and performance of membrane reacfor-Hydrogen separation membranes were 
fabricated by tape casting and rolling end-capped tubular supports (Ni-SrCe0 3 or Ni- 
SrCeo.8Zr 0i2 0 3 ) and by slurry-coating to form dense SrCeo.9Eu 0 .i0 3 . 5 (10ESC) hydrogen 



Figure 1. Integration of proton transport 
membranes (including SEM of membrane) 
in conversion of hydrocarbon fuels to H 2 . 
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permeable membrane thin films on the inner-side of tubular-type supports. In this work, the 
SrZr 0 . 2 Ceo. 7 Euo.i 0 3 . 5 tubular membrane, with Ni-SrCeo. 8 Zr 0 . 20 3 _ s porous support, was also 
investigated for the chemical stability improvement. Steam reforming and WGS reaction were 
carried out with this membrane. The theoretical CO conversion and hydrogen production were 
calculated using thermodynamic data. The significant effect of the membrane reactor on the 
steam reforming and WGS reaction is discussed. 

Experimental 

Fabrication and Characterization 
Powder Synthesis and Characterization 

Polycrystalline SrCe^xEuxO^ (x = 0.1, 0.15, and 0.2) samples were prepared by solid-state 
reaction. SrC0 3 (99.99%, Alfa Aesar), Ce0 2 (99.99%, Alfa Aesar), and Eu 2 0 3 (99.99%, Alfa 
Aesar) in the desired stoichiometric ratio were ball-milled for at least 24 h and then calcined for 
10 h at 1573 K in air. The calcined oxides were pressed into pellets and cold-isostatic-pressed. 
The pellets were sintered at 1723 K for 5 h in air and X-ray diffraction was used to confirm a 
single phase with the orthorhombic perovskite structure. 

Conductivity 

Pt electrodes were applied on the 10 mm diameter x 2 mm thickness dense disk pellets using 
Pt-paste (Engelhard 6926) and heated to 1273 K for 1 h for electrical measurement. 
Conductivity measurements were performed with a Solartron 1260 Impedance Analyzer in the 
frequency range of 0.1 Hz to 1 MHz and temperature range of 773 to 1 173 K under dry/wet pure 
hydrogen atmosphere. Water vapor was obtained by bubbling the hydrogen gas through de- 
ionized water at room temperature. 

Open circuit voltage (OCV) measurements were performed on dense disks 24 mm diameter x 
1.7 mm thickness. Each disk was sintered at 1723 K for 5 h with 2 °C/min ramp rate. The planar 
surface of each disk was connected with two platinum lead wires to measure OCV. The 
samples were sealed to the alumina tube with a glass sealant. 

The gas inlets in the set-up were positioned as close as possible to the test sample. A constant 
40 seem H 2 /Ar mixture was controlled by BOC Edwards 825 mass flow controllers on both gas 
inlets. Water vapor was obtained by controlling the bubbler temperature and gas tubes were 
wrapped using heating tape to prevent water vapor condensation. OCV was measured using a 
multimeter (Keithley 2000-20). 

Hydrogen Membrane Cell 

Ni-SrCe0 3 tubular-type supports for 10ESC hydrogen membrane thin film was prepared by tape 
casting as shown in Figure 2 and rolling techniques. Using tape casting and slurry coating, 
SrCe0 3 powder for a tubular-type support material was prepared by a conventional solid state 
reaction method from SrC0 3 (99.9%, Alfa-Aesar) and Ce0 2 powder (99.9%, Alfa- Aesar) as 
starting materials. The SrCe 0 . 9 Eu 0 .iO 3 powder was synthesized by a citrate process using 1 to 2 
molar ratios of the total metal nitrates to citric acid in order to prevent second phase formation 
and to decrease the calcining temperature [26], 
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Figure 2. Tape caster for making ceramic green tape. 
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Figure 3. Pictures of (a) tubular hydrogen membrane cells and (b) 10 mol% 
Eu-doped SrCe0 3 membrane coated on the inner side of NiO(Ni)-SrCe0 3 
tubular support at each processing step. 
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Figure 3(a) shows the 15 cm long end capped tubular-type hydrogen membrane cell. The 
preparation of the ESC hydrogen membrane on NiO (or Ni)-SC tubular support consists of four 
steps; and the pictures of a six-inch unit cell thus prepared at each preparation step are shown 
in Figure 3(b). 

Hydrogen Permeation Experimental Setup 

The SrCe 0 . 9 Euo. 1 C> 3 . 5 /Ni-SrCe 03 tubular-type membrane cell was installed in a hydrogen 
production reactor to confirm the cell was leak free. In order to make sure the leakage level and 
hydrogen permeation for the tubular-type hydrogen membrane cells, H 2 balanced with Ar was 
flowed to the feed side and He was flowed to the inner side of the cell to sweep permeated H 2 
and leaked Ar through the hydrogen membrane cell. The components of permeated gases on 
the sweep side were measured by the mass-spectrometer (Q100MS Dycor QuadLink). For the 
hydrogen permeation measurement, tubular-type hydrogen membrane cells were installed in a 
high temperature reactor apparatus as shown in Figure 4. One side of the membrane cell (feed) 
was exposed to H 2 (99.999%) diluted to the desired. 



Sweep gas (He) in 

Figure 4. Configuration of the hydrogen permeation 
reactor using a 15 cm long one-end-closed 
tubular-type hydrogen membrane cell. 
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Figure 5. Experimental set-up for hydrogen permeation test. 

concentration using He (99.999%) with 20~30 cm 3 /min total flow rate. For the wet gas flow, the 
feed gas was flow through a water bubbler at 25 °C to pick up 3 vol% water vapor. The other 
side (sweep) was swept with He at 20 cm 3 /min. Permeation was measured using a mass 
spectrometer (Dycor QuadLink IPS Quadrupole Gas Analyzer). The actual experimental set-up 
is shown in Figure 5. Permeated hydrogen was measured by a mass-spectrometer which was 
connected with the sweep gas outlet. Gas-chromatography was installed for measuring the 
product gases, while CH 4 gas was used as a feed gas. 

Results and Discussion 

Conductivity in SrCei. x Eux0 3 . 8 
Electrical Properties of SrCe^xEuxO^ss 

X-ray diffraction was used to determine the phase purity of SrCei. x Eu x 03 . s (x = 0.05 to 0.2). 
Single-phase with orthorhombic perovskite structure was achieved for each dopant 
concentration. The electrical conductivity of SrCei- x Eu x 0 3 (x = 0.1, 0.15, and 0.2), under various 
hydrogen atmospheres, was also investigated. In Figure 6, a plot of the electrical conductivity 
against P^ in dry hydrogen atmospheres shows a linear relationship between a tot and P ^ , 
which is consistent with Kosacki et al. for Yb-SrCe0 3 [21]. 

The hydrogen dependence of the conductivity can be explained if we assume that n-type 
conduction occurs because of the formation of electrons, as shown here 


H 2 (g) + 20 q <-> 20Hg + 2e' ; K, = = -^_ 

Hh 2 

(1) 

P H 0 [Vn ] n 2 

H 2 (g) + 0*^H 2 0(g) + V“ + 2e'; K 2 = H20t oJ 

K h 2 

(2) 
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(3) 


20Ho + 1 0 2 (g) o H 2 0(g) + 2h* + 20* ; 


RjjoP 

[OHo] 2 p | 2 



c 


c 

c 


Figure 6. Total conductivity of SrCe-i_xEu x 0 3 . 8 (x = 0.1, 0.15, and 0.2) as a function of P H 2 - 
The total conductivity of SrCeo.9Euo.-|0 3 _ 8 was higher than that of SrCe 0 .85Euo.i50 3 ^ and 
SrCeo.8Euo.20 3 . 8 , which are similar in value. Also, SrCe 0 .9Euo.i0 3 . 8 varied more with P H 2 than 
other formulations. This behavior may be attributed to protons being the predominant charge 
carrying species in SrCeo.9Eu 0 .i0 3 _ 8 , whereas SrCe 0 . 85 Eu 0 .i 5 O 3 . 8 and SrCeo. 8 Eu 0 .20 3 . 8 show non- 
negligible electronic conduction. 


Protonic, Electronic and Ambipolar Conductivity of SrCei. x Eu x 0 3 

Figures 7 and 8 show proton and electron conductivity as a function of dopant level (i.e., x in 
SrCe 1 . x Eu x 0 3 . 8 ) and temperature for strontium cerate in dry and wet H 2 , respectively. In both 
figures, proton conductivity decreases with increasing dopant concentration due to reduction of 
proton mobility and entropic destabilization of proton defects [29]. In contrast, electron 
conductivity increases with increasing Eu dopant. However, there is a slight decrease for x>0.15 
in SrCe-|. x Eu x 0 3 . 8 at low temperature due to the increased activation energy for electron 
conduction. Still, the decrease of electron conductivity in wet H 2 (Figure 8) was greater than that 
in dry H 2 (Figure 7) due to the higher Po 2 of wet H 2 . 

Generally, hydrogen permeation flux across the oxide membrane can be calculated using the 
Wagner equation (equation 4), assuming that bulk diffusion is the rate-limiting step. If oxygen 

conduction is negligible, then the proton flux, J QH ^ , is given by 
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Figure 7. Proton and electron conductivity behavior of SrCei-xEu x 0 3 -5 (x = 0.1, 0.15, and 0.2) 
with different temperatures under dry hydrogen atmosphere. 
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Figure 8. (a) Proton and (b) electron conductivity behavior of SrCe-i-xE^C^s (x = 0.1 , 0.15, and 
0.2) with different temperatures under hydrogen/water vapor atmosphere. 
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where a am b is the ambipolar conductivity, R is the gas constant, F is Faradays constant, a is 
conductivity and L is membrane thickness [20], 
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Equation 4 shows that hydrogen permeation is proportional to the ambipolar conductivity. 
Figure 9 shows the ambipolar conductivity dependence on dopant concentration in dry 
hydrogen atmosphere and hydrogen/water vapor atmosphere. The maximum ambipolar 
conductivity increases with temperature and Eu dopant concentration in dry hydrogen 
atmospheres. Therefore, the optimal dopant concentration for maximum ambipolar conductivity 
is obtained from 15 mol% Eu-doped strontium cerate at 600 °C to 20 mol% Eu-doped strontium 
cerate at 900 °C. Ambipolar conductivity values in hydrogen/water vapor atmosphere were 
lower than that in dry hydrogen atmosphere, and the maximum ambipolar conductivity was 
observed on 20 mol% Eu-doped strontium cerate in the range of investigated temperature. 
Therefore, maximum hydrogen flux can be achieved by controlling Po 2 , dopant concentration, 
and temperature. 



10 15 20 



mol % mol % 

Figure 9. The ambipolar conductivity behavior depends on the dopant concentration under 
(a) dry hydrogen and (b) hydrogen/water vapor atmospheres. 


Stability During Tube Fabrication of Hydrogen Membrane Cell 


Evaluation of Supporting Materials: NiO-SC 

Undoped strontium cerate was investigated as a support material. The main complication using 
SrCe0 3 as the support material is the reaction between SrCe0 3 and NiO at 1450 °C. Thus, NiO- 
SrCe0 3 tapes prepared at temperatures below 1450 °C were characterized. Observed 
shrinkage behavior of uniaxially-pressed SrCe^xEUxOs-s, SrCe0 3 and tape-casted SrCe0 3 -NiO 
shows that SrCei. x Eu x 0 3 . 5 sinters more easily than SrCe0 3 below 1450 °C and that SrCe-i. 
x Eu x 0 3 . 5 and NiO-SrCe0 3 have comparable shrinkage rates, which is necessary for dense thin 
film preparation. Moreover, the low sinterability of SrCe0 3 should improve the porosity of the 
support tube. Also, the samples with higher europium concentration demonstrate similar 
shrinkage rates. 

The thermal expansion of SrCei. x Eu x 0 3 -5 and 30 wt% NiO-SrCe0 3 must be similar in both air 
and in hydrogen to avoid failure of the SrCe 1 . x Eu x 0 3 . 5 film since the samples are sintered in air 
and tested in hydrogen. The expansion of SrCei. x Eu x 0 3 .5 is very similar to that of 30 wt% NiO- 
SrCe0 3 in air as shown by Figure 10(a) and the SrCe-|. x Eu x 0 3 -5 films do not crack after sintering 
which is confirmed by the SEM micrograph. In hydrogen, the thermal expansion coefficient 
(TEC) of the 30 wt% NiO-SrCe0 3 is approximately 13% higher than the TEC of SrCe 1 . x Eu x 0 3 . 6 . 
As shown by Figure 10(b), this expansion difference is not significant enough to result in 
cracking of the SrCei- x Eu x 0 3 -5 film. 
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(a) (b) 

Figure 10. Thermal expansion behavior of SrCei- x Eu x 0 3 -5 and NiO- SrCe0 3 in (a) air and (b) H 2 . 

Sintering Temperature of ESC film on NiO-SC Support 

The tubular-type hydrogen membrane cell is composed of a dense SrCe 0 . 9 Eu 0 . 1 O 3 . 6 film and a 
Ni-SrCe0 3 catalytic support. To preserve the integrity of the tube, the shrinkage rates between 
the SrCe 0 .9Euo.i0 3 . 5 film and the support need to be matched so that cracking does not occur on 
SrCeo.gEucuOs-g film. Figure 11 shows the linear shrinkage rates for NiO-SrCe0 3 support, which 
was prepared by tape casting and Eu 2 0 3 doped and un-doped SrCe0 3 which were prepared by 
uniaxial pressing according to the firing temperature. 
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Figure 11. Shrinkage 
rates for 1 0 mol% Eu- 
doped and undoped 
SrCe03 and SrCe03 
with 30 wt% NiO green 
tape. Microstructural 
images are shown for 
10 mol% Eu-doped 
SrCe03 thin films on 
NiO-SrCe03 supports. 
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<») <b) 

Figure 12. Cross-sectional view of (a) coated SrCeo.9Euo.i0 3 - 6 film on partially sintered 
NiO-SrCe0 3 and (b) dense SrCeo.9Euo.-i0 3 - 5 film on a Ni-SrCe0 3 support after sintering and 
reduction in H 2 atmosphere. 


By doping with Eu 2 0 3 , the SrCe0 3 system was easily densified, and the shrinkage rate for 
SrCe0 3 was increased. For the NiO-SrCe0 3 tubular-type support, the shrinkage rate of the 
support was higher than that of the SrCe0 3 thin films due to the high concentration of organic 
vehicles. In order to increase the compatibility of the linear shrinkage of SrCe 0 . 9 Euo.i0 3 . 6 film and 
the support, the NiO-SrCe0 3 support was initially partially sintered at 1100 to 1200 °C for 4 h. 
The SrCeo.9Eu 01 0 3 . 6 hydrogen membrane film which and the partially sintered NiO-SrCe0 3 
support were then sintered between 1350 to 1450 °C. The microstructures of SrCeo. 9 Eu 0 .i0 3 . 5 
film on NiO-SrCe0 3 supports according to sintering temperature are shown in Figure 1 1 . From 
the SEM results, the SrCeo.gEuo/iO^ film was fully densified above 1400 °C. Also, Figure 12 
shows the cross-sectional view of the as-coated SrCe 0 . 9 Euo.i 0 3 .5 layer on the inner side of the 
partially sintered NiO-SrCe0 3 tubular-type support and the fully densified SrCeo. 9 Euo.i0 3 -5, with 
apparent thickness of 30 |am, on the Ni-SrCe0 3 support after sintering at 1450 °C and reducing 
in a hydrogen atmosphere at 900 °C. The porous substrate was achieved by reduction of the 
NiO-SrCe0 3 at 900 °C. 


From our investigation, the final sintering temperature should be below 1450 °C to prevent 
reactions between NiO and strontium cerate, and above 1400 °C to obtain dense films of Eu- 
doped strontium cerate. In addition, to obtain compatible shrinkage rates between NiO- SrCe0 3 
and SrCe^xEUxO^s film during sintering, pre-sintering temperatures for the NiO-SrCe0 3 support 
should be in the range 1150 to 1200 °C. Cracks were observed in the NiO-SrCe0 3 support 
when it was pre-sintered at temperatures below 1150 °C; optimal pre-sintering and sintering 
temperature ranges are given in Figure 1 1 . 
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Hydrogen Permeation Properties at Dry/Wet Hydrogen Atmosphere 

The presence of water vapor reduces the concentration of oxygen vacancies and increases that 
of proton defects. However, hydrogen permeability in humid hydrogen atmosphere is lower due 
to lower electron concentration resulting from the higher P 0 2 of the humid atmosphere. 

Figure 13 shows the hydrogen permeation flux 
for H 2 (Ph2=0.03 atm) with and without water 
vapor in the feeding Ar gas (total 22 cc/min of 
flow rate) obtained from the hydrogen permeation 
reactor. As shown in this figure, hydrogen 
permeation in wet H 2 is slightly lower than that in 
dry H 2 . After reduction in H 2 , protons and 
electrons are the dominating defects. Under this 
condition, the hydrogen flux is expected to be 
proportional to [Ph 2 ] 1/4 according to Norby and 
Larring [22], 
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Figure 13. Hydrogen flux in dry/humid 
hydrogen as a function of temperature. 


The water vapor flux increased with increasing 
hydrogen partial pressure. In low oxygen partial 
pressures, oxides can be reduced, releasing 
oxygen which will form water vapor in the 
presence of H 2 . If the water vapor comes from 
the reaction of the permeated H 2 and oxygen 
ions in SrCe^xEUxC^ lattice, then the H 2 
permeation should be decreased as a result. 
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, , , , . Figure 14. Hydrogen flux as function of 

When the hydrogen permeation membrane is p d temoerature 

exposed to a wet hydrogen atmosphere, protons H2 
and electrons are the dominating defects and 
the hydrogen flux through the membrane is 
proportional to [Ph 2 ] 1/4 as shown in Figure 14. 
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Figure 15. Time dependence of Hydrogen flux according to the P H 2 with (a) Dry and (b) Wet H 2 . 


Figure 15 shows dry and wet hydrogen permeation flux on SrCeo.9Euo.i0 3 - Sl . The permeability of 
SrCeo. 9 Euo.i 0 3 . 8 under dry hydrogen atmosphere is higher than under wet hydrogen atmosphere. 
Under dry hydrogen atmosphere, hydrogen permeation flux through SrCeo.9Eu 0 .i0 3 - 8 degraded 
with time and this degradation was more severe under high hydrogen partial pressure on the 
feed side. On the contrary, the SrCe 0 .9Euo.i0 3 . 8 , under wet hydrogen atmosphere, shows stable 
hydrogen permeation flux with time. This degraded hydrogen permeation under dry hydrogen 
atmosphere may be interpreted as evidence for the secondary Ce0 2 phase formation by 
decomposition. 

Hydrogen Permeation Properties from CH 4 Steam Reforming 

We have evaluated the hydrogen permeability for the SrCei- x Eu x 0 3 .5 tubular-type hydrogen 
membrane cell in a H 2 atmosphere on the feed side. The overall goal of the SrCei- x Eu x 0 3 -5 
hydrogen membrane cell is pure hydrogen production from hydrocarbon fuels. Figure 16 shows 
hydrogen flux for CH 4 on the feed side. 

In order to use steam reforming of methane for hydrogen production using a hydrogen 
membrane cell, we have to consider methane conversion rates, partial pressure of product 
hydrogen, and solid state carbon formation rates. Using thermodynamic data for CH 4 and H 2 0, 
we calculated the conversion rate of methane, H 2 partial pressure after methane conversion, 
and the amount of carbon formation according to the ratio of methane to steam. When the 
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Hydrogen flux [cc/min] 


steam ratio is increased, the conversion rate of methane and the hydrogen partial pressure are 
increased. As shown in Figure 18, the hydrogen flux through the SrCei- x Eu x 0 3 -5 membrane for 
methane was higher than when pure H 2 was introduced. This is due to the increase of hydrogen 
partial pressure after methane was converted to CO, C0 2 , and H 2 by the steam reforming 
process. 


10 - 


18% humid methane 

(CH 4 8ccm and balanced Ar lSccm) 

3% humid hydrogen 

(H 2 8ccm and balanced Ar 15ccm) 


Nevertheless, the issue of coking needs to be considered when was using the hydrogen 
membrane cells on Ni-based porous support materials for the steam reforming of the methane 
process. Coking can cause degradation of hydrogen permeability and crack formation in the 

hydrogen permeation cells. As the steam 
ratio for methane increases, the formation of 
carbon decreases and the hydrogen 
permeation through the SrCe 1 . x Eux0 3 . 6 
membrane increases. For the steam 
reforming of methane, the methane to 
steam ratio should be below 0.6 to avoid 
coking and to increase the production of 
. .... pure hydrogen. 

b • — > 

Figure 17 shows the SrCe 0 . 9 Eu 0 .iO 3 . 6 /Ni- 
SrCe0 3 and SrCeo.gEuo.-iCWNi- 
SrCeo.8Zr 0i2 0 3 hydrogen membrane cells 
after exposed in methane with 8% steam. In 
the case of SrCe 0 . 9 Eu 0 .iO 3 _ 6 /Ni-SrCeO 3 cell, 
the crack occurred at around the 700 °C 
region. However, there were no cracks on 
the SrCe 0 . 9 Euo.i 03 - 5 /Ni-SrCeo. 8 Zro. 2 0 3 

hydrogen membrane cell. The crack is 
possibly caused by the phase change of 
SrCe0 3 in C0/C0 2 atmosphere and/or the 
coking in methane steam reforming that is above 0.6 ratio of CH 4 /H 2 0. Substitution of Zr onto 
Ce-site enhances the structural stability of SrCe0 3 and hence decreases the probability of 
cracking for the hydrogen membrane cell. Thus, the fabrication process for the tubular-type 
hydrogen membrane cells was modified accordingly. 
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Figure 16. Temperature dependence of 
H 2 flux in humid hydrogen/methane. 
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Figure 17. Photograph of the SrCeo. 9 Euo.i 0 3 - 5 /Ni-SrCe 0 3 and 10ESC/Ni-SrCe 0 .8Zr 0 .2O 3 hydrogen 
membrane cells after exposed in methane with 8% steam. 


IQESC/Ni-SrCeOa/Ni-SrCenRZrOnsOs Tubular-type Hydrogen Membrane Cells 

In order to prevent the reaction due to Zr-ion diffusion to SrCeo.9Euo.i0 3 - 5 layer, the Ni-SrCe0 3 
buffer layer was placed between SrCeo.9Euo.i0 3 - 6 and Ni-SrCeo. 8 ZrOo. 2 O 3 layer as shown in 
Figure 18. However, the Ni-SrCe0 3 buffer layer looks dense after sintering the hydrogen 
membrane cell. The hydrogen flux for SrCeo.9Eu 0 .i0 3 - 5 /Ni-SrCe0 3 /Ni-SrCeo.8ZrOo.20 3 hydrogen 
membrane cell was higher than that of the cell without the Ni-SrCe0 3 buffer layer, Figure 19, but 
the flux level was still lower than that of the SrCeo.gEuo.iCWNi-SrCeOs cell. To increase the 
hydrogen flux through the SrCeo. 9 Eu 0 .i 0 3 - 5 /Ni-SrCeo. 8 ZrOo. 20 3 hydrogen membrane cell, it is 
necessary to increase the porosity of Ni-IOESC buffer layer, thereby increasing the gas flow to 
the SrCeo. 9 Eu 01 03. 5 layer. 
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(a) (b) 

Figure 18. (a) FESEM image and (b) photograph of 10ESC/Ni-SrCe0 3 /Ni-SrCeo.8ZrOo. 2 0 3 
hydrogen membrane cell. 


Ni-SrCen«ZrQn?03 Tubular-type Supports 

To determine the phase stability of the cell supports in C0/C0 2 produced by methane steam 
reforming, the Ni-SrCeo. 8 ZrOo. 2 O 3 system was used instead of the previous Ni-SrCe0 3 systems. 
Figure 19 shows the hydrogen flux through the SrCeo. 9 Eu 01 03. 5 /Ni-SrCe 0 .8Zr0 0 . 2 O 3 hydrogen 
membrane cell. However, the flux level is much lower than that of SrCe 0 . 9 Euo.i 0 3 . 5 /Ni-SrCe 0 3 
cell. The SEM images show the grains of SrCe 0 . 9 Euo.i 0 3 - 5 close to the interface of 
SrCeo.gEuojOs-s and Ni-SrCe 0 .8Zr0 0 . 2 O 3 are larger than those close to the top of SrCe 0 . 9 Euo.i 0 3 _ 5 
layer. We supposed that the Zr ions reacted at the interface between the support and 
SrCe 0 . 9 Eu 0 .iO 3 .5 layer. 



P 1/4 [atm 1/4 ] 

H2 1 J 


Figure 19. Comparison of hydrogen 
flux for the SrCeo.gEuo.iCWNi- 
SrCeo. 8 ZrOo. 2 O 3 H 2 membrane 
cells with or without Ni-SrCe0 3 
buffer layer as a function of P H 2 1/4 
in humid hydrogen at 900 °C. 
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10ESZC or 15ESZC / Ni-SrCe n «Zr0n?O^ Tubular-type Hydrogen Membrane Cells 

Inevitable Zr ion diffusion into the membrane material resulted on Zr substituted ESC membrane 
materials. Though permeation flux could decrease by Zr substitution, continuous stable 
hydrogen permeation flux and chemical stability under C 0 /C 0 2 condition can be achieved. 
SrZr 0 .2Ceo.7Euo.i0 3 -5 and SrZro.2Ce 0 .65Euo.i50 3 -5 hydrogen membranes show increased hydrogen 
flux with temperature and hydrogen partial pressure. They also show [Ph2 ] 1/4 dependence of 
hydrogen permeation as in Norby and Larring’s model. Both hydrogen permeation show linear 
dependence with membrane thickness which means bulk diffusion is the rate-limiting factor for 
permeation. 

Water vapor can be formed from the reaction with permeated hydrogen gas and oxygen gas 
which comes from oxide lattice under reduced atmosphere. This water vapor formation may 
reduce the permeated hydrogen pressure. The permeated hydrogen from the both proton 
membrane (SrZro.2Ceo.7Euo/iO3_5 & SrZr 0 .2Ceo.65Eu 0 .i50 3 . 6 ) with counting on the hydrogen from 
water formation on the permeated side are shown in Figures 20 (a) and 20 (b), respectively. 
Increased hydrogen permeation with Eu dopant concentration in membrane material was 
observed. 

Figure 21 shows the hydrogen flux versus time at 5 % CO and 3 % FI2O at 900 °C. It decreased 
with time and an 8% drop was observed in 8 h for SrCe 0 .9Euo.i0 3 . 5 membrane. The hydrogen 
flux was stable for the SrZr 0 . 2 Ceo.7Euo.i0 3 -5 membrane, which verifies the stability enhancement 
by zirconium dopant. 




1/thickness 


Figure 20 . Permeated hydrogen from A) SrZr 0 .2Ceo.7Euo.i0 3 . 6 and B) SrZro. 2 Ceo.65Euo.i503- 5 with 
counting hydrogen from water formation. 
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Figure 21 . Hydrogen flux as a function 
of time under WGS reaction 
atmosphere. 
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SrCe0 3 -based Hydrogen Transport underwater Gas Shift Reaction 

After the supported hydrogen membrane was synthesized, a WGS reactor was fabricated by 
installing the hydrogen membrane to a quartz reactor and connecting it to the gas line by O- 
rings. The WGS reactor was put into a furnace vertically. The reactants (CO and steam) were 
flown from the top of the reactor to the outside of the hydrogen membrane. The retained gases 
were analyzed by gas chromatography (Varian CP 4900). A sweeping gas was flown to the 
inside of the hydrogen membrane. The outlet of sweeping gas was connected to a mass 
spectrometer so that the permeated hydrogen can be analyzed by the mass spectrometer 
(Q100MS Dycor, Quadlink). Gas flow rates were controlled by mass flow controllers. To 
investigate the effect of the WGS reactor on the hydrogen production and separation, three 
experiments were carried out with constant flow rates and concentrations of the reactants: 

(1) WGS reaction without hydrogen membrane: CO and H 2 0 were fed into the quartz 
reactor without installing the hydrogen membrane. 

(2) WGS reaction with hydrogen membrane, functioning solely as a catalyst: CO and H 2 0 
were fed to the outside of the hydrogen membrane with the sweep side of the 
membrane blocked to prevent hydrogen removal by permeation. 

(3) WGS reaction with the hydrogen membrane and in-situ hydrogen removal: CO and 
H 2 0 were fed to the outside of the hydrogen membrane. A sweeping gas was flown to 
the inner side of the membrane and the permeated hydrogen was carried out by the 
sweeping gas and analyzed by the mass spectrometer. 


The CO conversion and hydrogen selectivity is calculated as: 


X co % = 


pout 

— ^ x100% 


pout . pout 
°CO + °C0 2 


(5) 


X H2 % = - 


F s x Cl 


h 2 


: xC“'+F s xC ^ 2 


- x 1 00% 


( 6 ) 


NASA/CR— 2008-21 5440/PART3 


21 


where C^C^and C° u 2 ‘ are the concentrations of CO, C0 2 , and H 2 at the reactant outlet 
analyzed by the gas chromatography, respectively. C^is the concentration of permeated H 2 

analyzed by the mass spectrometer. F s and F R are the total flow rates of the sweeping gas and 
reactant gas outlet. In all three modes, two different gas compositions were applied: 3% CO 
+3% H 2 0 and 3% CO +6% H 2 0, while maintaining a constant flow rate balanced by Ar. Gas 
flow rates were controlled by mass flow controllers. 
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Figure 22. CO conversion versus temperature with (a) CO/H 2 0=1:1 and (b) C0/H 2 0=1:2. 
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Figure 22 shows the temperature dependence on CO conversion at three conditions with 
C0/H 2 0 ratio of 1:1 and 1:2, respectively. For the WGS without hydrogen membrane, CO 
conversion was lower than the thermodynamically calculated value, but increased with rising 
temperature for both C0/H 2 0=1 :1 and 1 :2. The reaction rate may be quite slow without catalyst. 
The dwell time for CO and H 2 0 in the reactor was very short. Equilibrium could not be reached 
before the gases exited the reactor, so the CO conversion was lower than the equilibrium value. 
On the other hand, the reaction rate increased with increasing temperature. Hence, a higher 
temperature resulted in higher reaction rate and higher CO conversion. With the catalyst, the 
CO conversion was slightly higher than the calculation value and decreased with increasing 
temperature. 


The WGS reaction consisted of a series of reactions, and the thermodynamic calculation was 
based on the assumption that all those reactions were in equilibrium. However, equilibrium may 
not have been achieved in such a short dwell time, even in the presence of Ni catalyst, which 
affected the CO conversion. Due to the exothermic nature of WGS reaction, CO conversion was 

more favorable at lower temperature. 
Figure 22(b) shows the CO 
conversion increased with 
increasing temperature with 
hydrogen in-situ removal and was 
higher than that with the catalyst. 
Timely in-situ removal of hydrogen, 
through the hydrogen membrane, 
reduced the hydrogen concentration 
on the product side and moved the 
reaction forward toward product side, 
resulting in a higher CO conversion. 

The hydrogen permeability of 
SrCeo. 9 Euo.i 0 3 . 5 increases with 
increasing temperature due to 
higher protonic and electronic 
conductivity [30]. A larger amount of 
H 2 permeated through 

SrCeo. 9 Euo.i 0 3 . 5 at higher 
temperature, leading to a larger 
improvement of CO conversion 
compared to that with the catalyst. 
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Figure 23. H2 flux versus temperature with a total flow 


rate of 20 ccm. 

Figure 23 shows, for the 

thermodynamic calculation, that 

hydrogen flux decreases with increasing temperature, in accordance with the CO conversion in 
Figure 22(a). The plot also shows the hydrogen production for WGS in the absence of a 
hydrogen membrane, increases with increasing temperature, a similar trend with the CO 
conversion, which can be explained by the kinetic reaction rate. However, without hydrogen 
removal hydrogen production decreased with increasing temperature, in agreement with 
(though slightly lower than) the trend of the calculated data. The difference was possibly caused 
by the by-reactions, such as the formation of methane, consuming some of the produced 
hydrogen. 
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Finally, Figure 23 shows the hydrogen permeation flux increased with temperature due to higher 
ambipolar conductivity at higher temperature. Likewise, the total hydrogen production increased 
with increasing temperature, in agreement with the trend of the CO conversion. At 900 °C, there 
was 39% improvement in hydrogen production compared to the equilibrium value. 

Conclusion 

The total conductivity of SrCe^EUxOs-s (x = 0.1 to 0.2) was measured by impedance 
spectroscopy, and ambipolar conductivity was calculated with proton and electron transference 
numbers. In the temperature range studied, SrCeo.9Eu 0 .i0 3 - 8 showed the highest total 
conductivity in both conditions. However, the proton transference number decreases with 
increasing europium dopant concentration, and electron transference number increases with 
increasing dopant concentration due to increasing n-type electronic conduction. In other words, 
proton conductivity decreases with increasing dopant concentration, while electron conductivity 
increases. A slight decrease in electronic conductivity for x>0.15, at low temperature, is 
observed due to the increased activation energy for electron conduction. The highest ambipolar 
conductivity was obtained from SrCe 0 .85Euo.i50 3 . 8 to SrCeo.8Eu 0 .20 3 . 8 between 600 °C and 
900 °C under dry hydrogen atmosphere and SrCe 0 .8Euo.20 3 .s under hydrogen/water vapor 
atmosphere. 

NiO-SrCe0 3 tubular-type supports for thin film hydrogen separation membranes were prepared 
by tape casting and rolling after the optimized slurry for tape casting was achieved. 
SrCe 0 . 9 Euo.i 0 3 . 8 thin film hydrogen separation membranes on the inner side of the partially 
sintered NiO-SrCe0 3 tubular supports were obtained by the colloidal coating method. Leak-free, 
15 cm long, end-capped, tubular SrCeo.gEuo.iOsyNi-SrCeOs hydrogen production cells were 
prepared by adjusting pre-sintering and final sintering temperatures. Hydrogen separation using 
the SrCeo. 9 Eu 0 .i0 3 . 8 /Ni-SrCe0 3 hydrogen membrane cell was successfully performed without 
leakage. 

Hydrogen production was improved significantly by incorporating the SrCe 0 .9Euo.i0 3 . 5 /Ni-SrCe0 3 
hydrogen separation cell for both the water-gas shift reaction and steam reforming. For the 
WGS reaction an improvement of 39% in hydrogen production compared to the equilibrium 
value was achieved at 900 °C. 
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Abstract 

Gas permeability in composite laminates is investigated based on Darcy’s law for porous 
materials. The permeability is derived in terms of transverse crack densities in each ply, 
delamination intersection area and a material constant to be determined from experiments. The 
permeability in cross-ply laminates is investigated first. A three-dimensional finite element model 
based on strain energy release rate is used to determine the crack densities in each ply of the 
laminate. Then the three-dimensional model with delaminations is used to obtain the crack 
intersection area. Parametric studies were performed and it has been found that the 
permeability is related not only to the resultant forces and thermal loads but also to length and 
shape of the delamination, stacking sequence and temperature dependent material properties. 
Finally, the progressive permeability in cross-ply laminate is predicted as a function of applied 
loads. 

Permeability tests were performed on various composite laminates including textile composites. 
The gas permeability, in general, increased with cryo-cycling. The increase was significantly 
higher in thin laminates. Textile composites showed remarkable damage tolerance in that the 
permeability remained al most constant after several cryo-cycles. The use of textile composites 
such as woven composites for cryogenic storage tanks warrants further investigation. 

Background 

High strength/weight and stiffness/weight ratio makes graphite/epoxy composite materials much 
more favorable to other conventional materials in aerospace industry. Reduction of structural 
weight of space vehicle is a critical issue to reduce the cost of sending payload to space. It is 
estimated that various gas storage tanks account for about 50% of the dry weight of space 
vehicles. Currently, the cost to transport payload into orbit is $10,000 per pound. NASA initiated 
the next generation single-stage-to-orbit reusable launch vehicle (RLV) program trying to reduce 
the cost to $1 ,000 per pound. One of such vehicle is the Lockheed Martin X-33 (see Fig. 1 ). 

The linear aerospike engines in X-33 operate using hydrogen mixed with an oxidizer, which 
needs to be stored at cryogenic temperatures. It had been estimated that composite fuel tank 
would reduce its weight by 40% and the overall weight of the vehicle by 14% compared to 
conventional metallic fuel tank (Kessler, Matuszeski, and McManus, 2001). The honeycomb 
sandwich structure was used in liquid hydrogen (LH2) tank in X-33. In November 1999, the 
outer wall of the sandwich structure explosively debonded when the tank was ground-tested at 
the NASA Marshall Space Flight Center. The reason for the failure was believed to be the 
excessive permeation of liquid hydrogen from inner face sheet to honeycomb. When the tank 
was warmed up, the gas in honeycomb could not escape quickly, and the corresponding 
pressure increase caused the failure. The failure has caused a setback to further application of 
composite materials in cryogenic tanks. Hence, understanding the mechanism of gas leakage 
through composite laminates is very critical to the future application of composite materials in 
cryogenic storage tanks. 

Composite laminate usually include several layers with different fiber orientations. Generally, the 
cryogenic tank is under high mechanical/thermal loading and also subjected to cyclic loads, 
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which could cause microcracks and delaminations. From experimental observation (Berthelot, 
2003), microcracks developed in 90°-ply of cross-ply laminate when the laminate was loaded in 
0°-direction. The “microcracks” in the present work refers to those matrix cracks fully developed 
in thickness direction of the ply and fully extended through the laminate. These microcracks 
induce local stress concentrations at the crack tip and can cause delaminations between 0°- 
and 90°-plies (see Fig. 2). 



Figure 1 : Lockheed Martin X-33 

(http://www.makelengineering.com/dir/Company/History/History.htm) 

In reality, the composite laminate is under bi-axial cyclic loading as the wall of the cryogen tank. 
In this case, the matrix microcracks will form in both 90°-plies and 0°-plies of cross-ply 
laminates. Because of the interlaminar delamination, the microcracks in different plies can be 
connected by intersection areas (shaded areas) as shown in Fig. 3. The network formed by 
microcracks and intersection areas provide a pathway for hydrogen to leak through the 
laminate. Of course, there will be some gas leak due to diffusion mechanism even no 
microcracks are present. But the amount would be negligibly small compared to the amount 
leaked through the cracks and delamination. Hence, in the present work, we focus on the gas 
leakage due to the damage - transverse cracks and delaminations - in the composite laminate. 




Figure 2: Microcracks in 90° ply and interfacial delaminations between 0°- and 90°- plies 
generated under thermal/mechanical loading. 
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Figure 3: Gas permeation pathway provided by microcracks and delaminations 
in cross-ply laminates. 

Literature Review 

There are a number of experimental and analytical investigations on microcracking formation in 
the literature, e.g., see Nairn (2000) and Berthelot (2003) and references therein. In the current 
research, we focus on the analytical modeling of microcracking. In conventional stress analysis 
it is assumed that microcracks form when stress in the 90°-plies reaches the strength of matrix 
material. This method, however, did a poor job to predict the microcracking process compared 
to the experimental observation (Nairn, 2000). Nairn presented methods based on fracture 
mechanics and energy balance to analyze microcracking, which assumed that new microcracks 
formed if the total energy released due to the formation of the microcracks reached the 
microcracking fracture toughness G mc . These energy release rate based methods did a good job 
of predicting microcracking in cross-ply laminates under uniaxial loading. The energy release 
rate method was then extended by Bapanapalli, Sankar and Primas (2006) to predict crack 
density evolution in cross-ply laminates under bi-axial and thermal loading using three- 
dimensional finite element method. This method is used in the current research to predict the 
crack densities evolution in each ply of composite laminates. 

After the failure of the liquid hydrogen tank of X-33 reusable launch vehicle many experimental 
and analytical works have been done to investigate the permeability in composite laminates. 
Stokes (2002) performed permeability testing on IM7/BMI laminated composites under bi-axial 
strains based on ASTM standard (ASTM D1 434-82, 1992). The permeability was measured 
under constant strain at room temperature and it has been found that the permeability is a time 
dependent property. Nettles (2001) tested permeability of composite laminates after impact 
testing and found that the flow rate usually had a nonlinear dependence on the applied 
pressure. McManus, Faust, and Uebelhart (2001) investigated the influence of loading 
conditions, crack density and ply orientation on the permeability of graphite/epoxy laminates. 
Choi (2005) performed permeability tests on various composite material systems subjected to 
certain cryogenic cycles, and found that for the laminates with same thickness the one with 
grouped lay-up has higher permeability than the one with dispersed lay-up. Among all the 
specimens, textile composite yields the lowest permeability. Kumazawa, Aoki, and Susuki 
(2003) investigated Flelium gas leakage through composite laminates, and compared the 
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experimental results with the analytical models. In their later work (Kumazawa, Susuki, & Aoki, 
2006), they measured the gas leakage rate of cross-ply laminates under biaxial loading and 
found that besides the amount of damage, the gas permeability is dependent on both load level 
and biaxial load ratio. Grenoble and Gates (2006) performed permeability tests on composite 
material IM-7/977-2 after specimen being mechanically cycled at room temperature and found 
that the leak rate depended on applied mechanical strain, crack density and test temperature. 
Bechel, Camping, and Kim (2005) measured the evolution of microcrack density in composites 
due to cryogenic thermal cycling. Bechel and Arnold (2006) investigated the permeability in 
multi-directional laminates after certain thermal cycling from room temperature to cryogenic 
temperature and found that the permeability increases significantly as the number of cycling 
increased. 

Besides the experimental measurement of the permeability in composite laminates, efforts were 
also devoted to develop analytical and numerical models to understand the mechanism of the 
permeability. As shown in Fig. 3, the matrix microcracking and the interlaminar delaminations 
form the leakage pathways for cryogenic gas. The intersection areas (the shaded areas shown 
in Fig. 3) formed by the crack opening in adjacent plies are the neck areas in the leakage 
pathway. If we treat the damaged composite laminate as a porous medium, then Darcy’s law for 
porous materials can be used in modeling gas permeability (Aoki, Kumazawa, & Susuki, 2002, 
Kumazawa, et al. 2003, Roy & Benjamin, 2004, and Roy & Nair, 2006). In Darcy’s law, the 
volume flow rate of fluid is proportional to the pressure gradient, and the proportionality constant 
is given by the ratio of permeability constant to the viscosity of the fluid. The gas permeability 
constant is determined by the micro-structure of the porous material. In our case, it is 
determined by crack densities in each ply, ply thicknesses, intersection area and a material 
constant that can be determined from experiments. For a given laminate, the ply thicknesses 
and the material constant are fixed. The crack densities can be predicted using the energy 
release rate method mentioned before or directly measured from experiments. However, the 
crack opening of the microcracks is hard to measure in experiments. Hence, analytical models 
are required to predict the intersection area. Kumazawa et al. (2003) only considered 
microcracking in their model and took mean crack opening displacement computed from two- 
dimensional shear-lag analysis to get the intersection area. Theoretically, however, 
microcracking alone cannot form the intersection area; delaminations have to be introduced in 
the model. Roy and Benjamin (2003) developed a model with delamination and computed crack 
opening displacement based on first-order shear deformable laminate theory. They compared 
the analytical results for crack opening displacement with finite element results and investigated 
the permeability under mechanical loading and thermal loading. Later on, Roy and Nair (2006) 
extended this model to include the stitch cracks developed in the oblique layers of the laminate 
by introducing effective spring stiffness. All these models are 2-D models. In the current work, 
we have found that the intersection areas are largely dependent on the local geometries and a 
full three-dimensional analysis is required to obtain accurate result (Xu & Sankar, 2007). 

Some researchers predicted the gas permeability in composite laminate based on fluid 
mechanics principles (Noh, Whitcomb, Peddiraju, & Lagoudas, 2004, Peddiraju, Lagoudas, 
Noh, & Whitcomb, 2004, Peddiraju, Grenoble, Fried, Gates, & Lagoudas, 2005, and Peddiraju, 
Noh, Whitcomb, Lagoudas, 2007). In their models, the opening due to transverse microcracks 
and delaminations and the intersection area are first investigated. Then a damage network 
through laminates with assumed crack opening and delamination was set up and the effective 
conductance through this network was predicted using a computational fluid dynamics program 
FLUENT®. In these studies, the effective conductance is related to the ply thickness, the crack 
density in each ply. Also, the conductance of gas through a representative volume element 
(RVE) of intersection and microcracks were computed and a simplified model was set up, which 


NASA/CR— 2008-21 5440/PART3 


32 



agreed well with the previous simulation results. In these models, however, since the damage 
network was assumed with fixed opening and intersection area, each model can only represent 
one single loading case because the opening and the intersection area are dependent on the 
applied loads. Hence these models have difficulties to relate conductance to applied loads. 

Most of the analytical models focus on cross-ply laminates, in which transverse matrix cracks 
are developed in each ply and extended thorough the thickness in both thickness direction and 
width direction. In general, however, it is good to include angle-plies in composite laminate 
design to gain better performance. Different from extensive transverse matrix cracks in cross-ply 
laminates, short cracks (or called stitch cracks) emanated in angle plies in multi-directional 
laminates. This phenomenon was first identified by Jamison, Schulte, Reifsnider, and 
Stinchcomb (1984) in the multi-directional laminates under fatigue loading. Lately, the stitch 
cracks formed in laminates with stacking sequence [+0 n /-0 n /9O 2 n]s were systematically 
investigated by Lavoie and Adolfsson (2001). It has been found that the stitch cracks not only 
emanated under fatigue loading but also emanated due to monotonic loading and thermal 
residue stresses. Stitch cracks were also observed in [0/+45/-45/90] s laminates by Bechel et al. 
(2005) after certain thermal cycling between room temperature and cryogenic temperature. 
Yokozeki and his colleagues (Yokozeki, Aoki, Ogasawara, and Ishikawa (2005), Yokozeki & 
Ishikawa (2005), and Yokozeki, Aoki, & Ishikawa (2005)) performed experiments on laminates 
of type [O/0 n /9O] s under uniaxial loading and they found that whether stitch cracks or developed 
cracks will emanate in 0°-ply depends on the value of 9 and the thickness of the 0°-ply. 

Some of the analytical models mentioned before included stitch cracks in the modeling. As 
mentioned before, Roy and Nair (2006) introduced an effective spring stiffness to incorporate 
the stitch crack into their model and still used the first-order shear deformable laminate theory to 
compute the crack opening displacement of stitch crack. In the damage network, Peddiraju et al. 
(2005) included the stitch cracks with assumed crack length and crack opening. Note that all 
these models either eliminated the difference between stitch crack and transverse crack using 
an effective spring stiffness or assumed damage state of stitch cracks. They were not able to 
represent the relationship between stitch crack and loads. In reality, the crack densities, the 
crack length and the opening displacement of stitch cracks are functions of the applied loads. 

Objectives 

The objectives of the current work are: 

• Develop a method based on Darcy’s law and three-dimensional finite element analysis to 
predict the gas permeability in cross-ply composite laminates. 

• Extend the method to multi-directional laminates 

• Study the effects of ply thickness, delamination length, delamination shape, lay-ups on the 
permeability 

• Perform permeability tests on a variety of composite laminates including textile composites 

• Study the effects of cryo-cycling on the permeability 

Analytical Model 

There are several theories that might be useful in deriving the permeation model for composite 
laminates: Darcy’s Law for porous materials, Poiseuille’s Law for capillary and Darcy-Weisbach 
Equation for pipe pressure loss. In the present work, the permeation model in composite 
laminate is mainly derived from Darcy’s law for porous materials and verified by other laws. The 
permeability is represented in terms of crack density, crack intersection area and a material 
constant to be determined from experiments. The derivation of the permeation model is based 


NASA/CR— 2008-21 5440/PART3 


33 



on the cross-ply laminate, where only 0°-ply and 90°-ply present, but it can be extended to multi- 
directional-ply laminate as well. 

Permeation Model 

From experimental observations, matrix transverse micro-cracking in 90°-ply is the first form of 
damage in cross-ply composite laminates. Due to stress singularity at the crack tip when 
transverse crack reaches the interface between 90°-ply and 0°-ply, interlaminar delamination 
may occur at the interface. Theoretically, microcracks alone cannot form a contiguous path for 
gas leakage since no intersection area is formed. Hence the delamination has to be introduced 
in the numerical model as shown in Fig. 3. The shaded intersection areas formed by the crack 
opening displacements (CODs) of adjacent plies connect the micro-cracks, and allow gas to 
leak out. If the composite laminate with micro-cracks and delaminations were treated as a 
porous material, Darcy’s Law for porous materials can be applied. 

Darcy’s law for viscous flow of gases through porous media is given as 


q = -—VP 
A 


( 1 ) 


where q is volume flow rate per unite area, k is permeability tensor, p is the viscosity and VP is 
pressure gradient vector. If only consider the gas leak in thickness direction of laminate, 
following one-dimensional Darcy’s law can be used: 


k_dP_ 
/u dz 


( 2 ) 


where k is overall permeability constant in thickness direction and z is the thickness coordinate. 
The permeability constant k is dependent on the microstructure of the laminate. Integrating both 
side of Eq. (2) over laminate thickness gives 

qh = --(p N -P 0 ) (3) 


or 


-(P»-Po) = <7hf 


( 4 ) 


where h is the laminate thickness, P N and P 0 are pressure at the top and bottom respectively 
(see Fig. 4). 
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Figure 4: Pressure distribution along thickness direction of laminate. 

Also, if the integration is done layer-wisely, we have 

-(P,-Pm) = <7'i,p / = 1 N (5) 


where h\ and k\ are thickness and permeability for each lamina. Summing Eq. (5) from / = 1 to 
/ = N yields 


-(P N -P 0 ) = qMf,T 

i = 1 K i 


Compare Eq. (6) with Eq. (4), it can be seen that 

h A hj 
— = > — or 

k fak, 


k = 


h 


I 

/=l 


( 6 ) 


(7) 


For each lamina, most of the resistance comes from the intersection area. According to the 
Poiseuille’s Law, the volume flow rate in a rectangular tube (see Fig. 5) for the laminar flow is 
given by: 


Q = 


d 4 

12//L 


AP 


Dividing both sides by the cross-section area of the tube, d 2 , gives: 


d 2 AP 
12// L 


( 8 ) 


(9) 


Comparing Eq. (9) with Eq. (2), it can be seen that the permeability constant k for the 
rectangular tube is proportional to the cross-section area. Based on this fact, we can assume 
that the permeability for each lamina is proportional to the total intersection area. This 
assumption is also confirmed by Peddiraju, et al. (2007) using FLUENT® simulation. Noting that 
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first and last laminas only have one interface with other laminas while others have two 
interfaces, so the individual permeability is assumed to take following form: 


/c 1 — C/L,/l 2 A (1 2 ) 

ki = c(/i j _ 1 /i l A {j _ rr) + /l ( /l (+1 /A ((( ^ij)/ 2, i = 2,...,N — 1 
/c N = CA N _ 1 A N A (N _ 1N ^ 


( 10 ) 


where /l,- is the crack density in / <h ply, ^(/, /+ i) is the individual intersection area between / h and 
/+1 th ply, and C is a material property to be determined by experiments. 



Substituting ki back into Eq. (7), the overall permeability of laminate can be written as: 


k = 


Ch 


V^Al.2) + L=2 2 V( /l /-l' l /A/-1,/) + ^/^/+iAm+1)) + ^w/Vi^nAw-V 


( 11 ) 


-1,N) 


It is obvious from Eq. (11) that the permeability is determined by the crack densities A , and the 
intersection areas A^ once the laminate lay-up is given and the material constant C is 
characterized by experiments. In another word, the overall permeability is determined by the 
total intersection area in each ply because that the crack densities give the number of 
intersection area in each ply. 

Preliminary Two-Dimensional FE Modeling 

Before using a three-dimensional finite element model, two-dimensional FE modeling is 
presented first to obtain some preliminary results. A three-layer cross-ply laminate is considered 
in this modeling. Figure 6 shows the representative volume element (RVE) for two different 
cases of damage scenarios that could occur in cross-ply composite laminates. In the first case, 
[90/0/90], the delamination emanates from the surface ply of the laminate whereas in the 
second case, [0/90/0], it emanates from the middle ply. Due to symmetry, only one-fourth of the 
RVE is considered in the simulation. Eight-node isoparametric elements are used in the 
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simulation and the orthotropic material properties used in the simulation are shown in Table 1. 
The simulation has been done using the commercial finite element package ABAQUS 6.5®. 

From Fig. 6, it can be seen that the RVE can be described by crack spacing (reciprocal of crack 
density), ply thickness and delamination length. The ply thickness h is taken as 0.33 mm in the 
current model. First, the variation of strain energy release rate (SERR) at the delamination tip 
with respect to the delamination length a is investigated. The crack density is taken as 1 cm' 1 . 
Even though the individual Mode I and Mode II strain-energy release rates do not exist due to 
the oscillatory characteristic of stresses and displacements near the crack tip for bi-material 
interface cracks (Sun & Jih, 1987), the total strain-energy release rate is well defined 
(Hutchinson & Suo, 1992). Under linear elastic assumption, path independent J-integral is 
identical to strain energy release rate. Since J-integral is insensitive to the FEM mesh, uniform 
mesh is accurate enough for this simulation. 
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Figure 6: Two-dimensional delamination model: (a) Delamination from surface plies [90/0/90]; 
(b) delamination from middle ply [0/90/0]. 


Table 1: Orthotropic Material Properties for the Composites. Moduli are in GPa 


Properties 

E i 

e 2 =e 3 

Gl2— G-13 

G 23 

Vl2=Vi3 

V23 

OL\ 

a 2 = a 3 

Value 

169 

8.62 

5.0 

1.22 

0.355 

0.41 

-3.384e" 9 /°C 

28.998e" 6 /°C 


Figure 7 shows the J-integrals with respect to the delamination length a. The delamination 
length a is normalized by the ply thickness h in Fig. 7. The RVE is under constant displacement 
loading. It can be seen that the J-integrals reach maximum values when delamination length is 
about the ply thickness, and then decrease slowly. This implies that the delamination propagate 
steadily when delamination length exceeds ply thickness under displacement control loading. 
When the delamination length is too small, less than the ply thickness, the boundary effects are 
obvious. 

From Fig. 7, it seems to be appropriate to take ply thickness as the delamination length since 
the delamination propagates steadily after its length exceeds ply thickness. Hence, in three- 
dimensional FE analysis, a length comparable to the ply thickness is taken to be the 
delamination length and the sensitivity of the intersection area to the delamination length is 
investigated. 
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Figure 7: J-integral variation with respect to delamination length a, which is normalized by 
dividing by ply thickness h\ (a) [90/0/90] laminate; (b) [0/90/0] laminate. 

Next, the crack opening displacements (COD) of microcracks are computed with a fixed 
delamination length and varied crack density. Even though both crack density and intersection 
area are functions of applied loads, fortunately, following results show that these two factors can 
be decoupled since crack density does not have much influence on the COD. Figure 8 shows 
the computed COD with respect to crack density. The CODs are expressed for unit force by 
normalizing with respect to the force resultant N x , which can be computed directly from finite 
element analysis by volume average method (force resultant equal to the summation of product 
of stress and element volume and then divided by the laminate thickness). It can be seen from 
these two figures that the CODs are almost constants for a given force resultant and 
delamination length. 
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Figure 8: (a) COD variation with respect to crack density for lay-up [90/0/90]; (b) COD variation 
with respect to crack density for lay-up [0/90/0], 
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From the above discussion, we can argue that the crack density does not influence the crack 
opening displacement much. Hence, the permeability prediction can be decoupled into two 
parts, first find crack densities for given force resultants and then compute intersection area 
under these loads. Then, the permeability for a given set of force resultants can be computed 
using Eq. (11). The crack densities in each ply can be either predicted using analytical methods 
(Nairn, 2000 and Bapanapalli, et al., 2006) or measured directly from experiments (Nairn, et al., 
1993 and Nairn, 2000). However, most of the available experimental data are for laminates 
under uniaxial loading. In reality, the composite laminate in cryogenic tank is under biaxial 
loading. Hence, the analytical method developed by Bapanapalli, et al. (2006) will be adopted in 
the current work to predict crack densities in cross-ply laminate under biaxial loading. The 
intersection area will be computed using three-dimensional FEM and a parametric investigation 
will be performed. 

Crack Density Prediction 

In this chapter, the method developed by Bapanapalli, et al. (2006) is briefly summarized, 
predicted crack densities under uniaxial loading is compared with experimental result, and then 
crack densities under biaxial loading are predicted in various cases. 

A composite laminate with configuration [0/90/0] is considered in the current work. From 
experimental observation (Nairn, 1993), two different microcracking scenarios could develop in 
the laminate when it is subjected to uniaxial loading as seen in Fig. 9. If the laminate is loaded in 
0°-direction, microcracks emanated in 90°-ply (middle ply) and next new microcrack usually 
formed in between two existing cracks (see Fig. 9(a)). On the other hand, if the laminate is 
loaded in 90°-direction, microcracks emanated in 0°-plies (surface plies). The asymmetric crack 
pattern shown in Fig. 9(b) has been observed by Nairn (2000). Hence, four new cracks formed 
to keep the same pattern. 

A finite fracture mechanics approach was used by Nairn (2000) to predict new microcrack 
formation. In this approach, the strain energy release rate due to the formation of new 
microcracks is expressed as a function of applied load, mechanical properties and crack 
density. Then, this strain energy release rate is equated to microcracking fracture toughness, 
G mc , to solve for the required load to form new microcracks. It has been shown that the 
microcracking fracture toughness G mc is a material property and it can be determined through 
experimental techniques (Nairn, 2000). The predicted results using this approach matched the 
experimental results quite well. However, this approach can only predict the crack density of 
composite laminate subjected to uniaxial loading. Bapanapalli, et al. (2006) extended this 
approach to predict crack densities in both 0°-ply and 90°-ply at same time when the laminate is 
subjected to biaxial loading. 

In the current work, the extended finite fracture mechanics approach (Bapanapalli, et al., 2006) 
is used to predict crack densities in each ply of composite laminate. Based on the two scenarios 
shown in Fig. 9, a three-dimensional unit cell is taken as shown in Fig. 10. In this model, crack 
surfaces are shown in shaded area. The dimension of the unit cell is determined by the crack 
densities of surface and middle plies. 

The strain energy of the unit cell for load control condition is given by 


111 — 

G = — — — -N T A N 
2 /l x A y 


( 12 ) 
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where X = P' A 12 

A A 

™21 ^*22 


is the inverse of laminate extensional stiffness matrix (or the 


compliance matrix), N= is force resultant matrix, and X x , X y are crack densities in x- and 

_ y _ 

y-direction. The formation of new microcracks will cause the degradation of the stiffness matrix 
of the laminate. Hence, the released strain energy due to the formation of new microcracks is 
given by 

G m = — — — —N t AA N (13) 


where AA = A 2 - A, is the difference of the compliance matrix before and after the formation of 
new microcracks. 


En, Co En, Cn ft' Cn n A 
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(a) (b) 

Figure 9: Two microcracking scenarios in composite laminate [0/90/0] under uniaxial loading, (a) 
Unit cell for microcracks in middle ply and new crack formation, (b) Unit cell for microcracks 
in surface ply and new cracks formation. (Courtesy: Bapanapalli, et al., 2006) 
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Figure 10: Three-dimensional unit cell with microcracks in both middle and surface plies. 


As the wall of cryogenic tank, the composite laminate is under biaxial loading, and the hoop 
stress and the longitudinal stress are proportional to each other. Hence, it is reasonable to 
assume that the force resultants in the laminate in the form: 

N x =aN y (14) 

where a is proportional factor, whose value equals to 2 or 0.5 for the cylindrical body and 1 for 
the spherical cap of a pressure vessel. Substituting Eq. (14) into Eq. (12), the released strain 
energy can be rewritten as 

G m =l~(a 2 AA„ + 2aAA, 2+ AA 22 )N 2 y (15) 

Z A y 

where A A u , A A u , and AA 22 are the entries of matrix A A . 

According to the finite fracture mechanics, the strain energy released due to new crack 
formation should equal to the energy required to create new crack surfaces. There are three 
cases could happen for the microcrack formation: 

• Case I: new microcracks form only in the surface plies. In this case, X x remains 
the same but A y is tripled. 

• Case II: new microcrack form only in the middle ply. In this case, JL X is doubled 
and / L y remains the same. 

• Case III: new microcracks form in both middle and surface plies. In this case, A x 
is doubled and JL y is tripled. 

The required surface energies for these three cases are given by 
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( 16 ) 


U,=4j-h„G mc 


1 

U/i = -^—h ao G mc 

A y 

1 1 

U lll = ~T h QO G mc +4 ~T h 0 G mc 

A y A x 


where h 0 and h 90 are the thickness of 0°-ply and 90°-ply respectively. 

Equating these surface energies to the released strain energy in Eq. (15), one can solve for the 
required force resultants for the three cases as follows 


N 

y(i) „2 


8h 0 A y G mc 


cc AA U + -I - A^22 


n „», - 


2h 90 A x G mc 


cc AA^ + ^ccAA^ AA22 


N 


^mC ( 2^90 ) 


y(iii) 


cc + r 2-(xAJK^2 +A A 


22 


(17) 

(18) 


(19) 


In reality only one case could happen, so the case that gives the smallest value of A/ y will 
happen. In another word, for a given unit cell (A x and A y fixed), if N y( n) is the least values among 
three computed force resultant, then next microcrack will form in surface ply. Note that the 

values of A A are different for these three cases. 

This method can also be extended to the displacement control condition. For details, please 
refer to Bapanapalli, etal. (2006). 

Predicted Analytical Results 

Uniaxial Loading — First, the crack density in composite laminate is predicted under uniaxial 
loading condition. To compare with the available experimental data, laminates of type [0/90 n /0] 
(n=1, 2, 4) are investigated. The material properties are the same as in Table 1. 

In Eq. (17) - Eq. (19), for each combination of (A x , Ay), the stiffness matrix [A] is different. Matrix 
[A] is computed using FE analysis, but if FE analysis is performed for every unit cell (given A x 
and A y ), there will be too much computation since 3-D FE analysis is very time consuming. 
Hence, in the current work, the [A] matrices for certain (A x , A y ) combinations are computed and 
then stiffness matrix is fitted as a complete cubic function of A x and A y based on the computed 
values using least the square approximation. 

For a given unit cell, the [A] matrix is computed based on constitutive relation: 
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Displacement boundary conditions are applied to the unit cell to make the strain status be s x = 1 
and e y = 0. The force resultants N x and N y are computed using volume average method. Then 
the first column in [A] matrix equals to the computed force resultants. Similarly, applying the 
boundary condition s x = 0 and e y = 1 will give the second column of [A] matrix. For each pair of 
(A x , A y ), two FE analysis are needed to determine [A] matrix. In our case, seven A x (ranging from 
0 to 9.5 cracks/cm) and eight A y (ranging from 0 to 12.6 cracks/cm) are taken. So totally 1 12 FE 
analysis are needed for each laminate to determine response surface of [A] matrix. It can be 
seen that this is a quite time consuming procedure. Figure 1 1 shows the response surfaces for 
the entries in [A] matrix of laminate [0/90/0]. The dot marks are data points. Other laminates 
also have similar response surfaces. 

A MATLAB program is developed to predict the progressive crack densities A x and A y of 
laminate subjected to biaxial loading as a function of N x or N y . In this program, initial values of 
crack densities and the ratio between N x and N y . For uniaxial loading, the ratio a is taken as 0. 

From experimental observation (Nairn, 2000) in [0/90 n /0] laminates, it has been found that 
initiation stress for microcracks in middle ply is lower for thicker middle ply and initial stress for 
microcracks in surface ply is higher for thicker middle ply. However, at high stress level, larger 
crack density was observed in the middle ply with thinner middle ply. The predicted crack 
densities of laminates [0/90„/0] (n = 1, 2, 4) under uniaxial loading are shown in Fig. 12 and 
Fig. 13. The microcracking fracture toughness G mc is taken as 240 J/m. In Fig. 12, the stress to 
initiate first microcrack increases as the middle ply thickness decreases. However, once the 
initiation stress is reached, the crack density increases faster in the laminate with thinner middle 
ply. For [0/90/0], even a negative slope is obtained, which means that the microcracks 
accumulate simultaneously in the laminate. As seen in Fig. 13, as middle ply thickness of 
laminate increases, the initiation stress increase as well. Hence, our prediction and 
experimental observation matches quite well qualitatively, which gives us confidence to use this 
approach to predict the crack densities of laminates under biaxial loading. 

Biaxial Loading — In reality, the laminate used in propellant tank is under bi-axial loading. The 
proportional constant a between force resultants N x and N y is taken as 0.5, 1.0, and 2.0 to 
represent different possible stress states in hydrogen tank. Figure 14 to Fig. 22 show the 
variations of crack densities in different laminates under different biaxial loading. For each 
laminate, crack densities A x and A y are plotted in same figure with respect to the average stress 
in y-direction (Ny divided by the laminate thickness). Figure 14 to Fig. 16 represent the cases 
when a = 0.5, Fig. 17 to Fig. 19 represent a = 1 .0, and Fig. 20 to Fig. 22 represent a = 2.0. From 
these figures, it can be seen that the initiation stresses for first microcrack in both x- and y- 
direction are lower than their corresponding values under uniaxial loading. With the same 
proportional constant a, the initiation stress for A x decrease dramatically as the middle ply 
thickness increase. On the contrary, the initiation stress for A y increases as the middle ply 
become thickness, but the increase step is not as large as that under the uniaxial loading. 
Loads in x-direction affect the crack densities in y-direction. For the same laminate, say [0/90/0], 
when a increases (<j x increases), the initiation stress for A y decreases. No microcracks are 
generated in [0/90/0] for a = 0.5 and 1 .0 when A y already meet the preset limit. These predicted 
crack densities will be used to predict the progressive gas permeability later on. 
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Figure 1 1 : Response surface for [A] matrix of laminate [0/90/0]. (a). A1 1 , (b). A12, (c). A22. 
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Figure 12: Crack density evolution in middle ply under loading in x-direction. 
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Figure 13: Crack density evolution in surface ply under loading in y-direction. 
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Figure 14: Crack density prediction for laminate [0/90/0]. Nx/Ny = 0.5. 
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Figure 15: Crack density prediction for laminate [0/902/0]. Nx/Ny = 0.5. 
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Figure 16: Crack density prediction for laminate [0/904/0]. Nx/Ny = 0.5 
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Figure 17: Crack density prediction for laminate [0/90/0], Nx/Ny = 1.0. 



Figure 18: Crack density prediction for laminate [0/902/0]. Nx/Ny = 1.0 
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Figure 19: Crack density prediction for laminate [0/904/0]. Nx/Ny = 1.0. 
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Figure 20: Crack density prediction for laminate [0/90/0]. Nx/Ny = 2.0. 
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Figure 21: Crack density prediction for laminate [0/902/0]. Nx/Ny = 2.0. 
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Figure 22: Crack density prediction for laminate [0/904/0], Nx/Ny = 2.0. 


Investigation of Intersection Area 

After predicting the crack density, we will focus on crack intersection area and its variation due 
to different factors. 

Model Description — Three-layer cross-ply composite laminate [0/90/0] is considered in the 
current model and following assumptions have been made for simplicity: the crack densities and 
the position of microcracks in two surface plies are identical; the microcracks in both surface 
plies and middle ply are evenly distributed. Once crack densities are given, the dimension of the 
unit cell is fixed as shown in Fig. 23. The upper figure in Fig. 23(a) represents the top view of 
the laminate with transverse matrix microcracks. The horizontal and vertical solid lines represent 
the microcracks in 0°- and 90°- plies respectively. Due to both in-plane symmetry and symmetry 
in thickness direction, only one-eighth of the unit cell is analyzed as shown in bottom figure in 
Fig. 23(a). Two free surfaces are the microcrack faces. Symmetric boundary conditions are 
applied to the surfaces opposite to the free surfaces and bottom surface and constant strains 
are applied as shown by the arrows. Delamination is introduced in this model since no 
contiguous gas-leak-pathway will be formed without delamination. The delaminations emanated 
from both surface ply and middle ply form a rectangle delamination front. In the current model, 
the crack densities in the 0°- and 90°-plies are taken as equal, A x = A y = 1 cm" 1 . The material 
properties are kept identical to those of the two-dimensional model as given in Table 1. Ply 
thickness is still taken as 0.33 mm. Twenty-node three-dimensional isoparametric solid 
elements and uniform mesh are used in the simulation. About 28,000 elements were used in the 
model. Figure 23(b) shows a typical deformation of the model under bi-axial loading. Note that 
the deformation is enlarged by a factor of 10. 
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Figure 23: Three-dimensional delamination model and typical deformation. 

Kumazawa et al. (2003) and Roy and Benjamin (2004) proposed to superpose the two- 
dimensional CODs obtained from analytical methods to get the intersection area. In fact, if we 
monitor the COD along the path shown by the arrow in Fig. 24(a), the COD first remains as a 
constant at the straight edge, and then increases remarkably when close to the rectangular 
crack front. The constant value along the straight edge is identical to the result predicted by two- 
dimensional analytical model. The COD near the rectangular delamination front, however, is 
almost 50% larger than the two-dimensional result as shown in Fig. 24(b). Therefore, the local 
effects cannot be neglected and COD near the rectangular front computed from three- 
dimensional analysis should be used to compute the intersection area. Three-dimensional 
analytical model is too complicated, hence all the CODs in the current work are computed using 
FE analysis. 
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Figure 24: COD along specific path shows the local effect. 
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Intersection Area Computation — Due to the linear relation between COD and load, the 
intersection area under any loading case (A/ x , A/ y , AT) can be obtained once the COD for the 
following three basic cases are computed: 
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# (From room temperature to cryogenic temperature). 

For any combination of (A/ x , A/ y , AT), the corresponding (£ x , £ y , AT) can be obtained from the 
basic formula for composite laminate shown in Eq. (21), in which the matrices [A] and [a] can be 
determined from FE analysis of those three basic cases using the approach mentioned before. 



The corresponding CODs for (Nx, Ny, AT) are then obtained by superposing the results of three 
typical cases with different weighting factors. The intersection area AA is computed by 
assuming that the cross-section formed by the CODs at corner is a rectangle. Figure 25 shows 
the variation of the intersection area with respect to the mechanical and thermal loads. In 
Fig. 25(a), different symbols represent different ratios between Ny and Nx, the solid lines 
represent the cases with thermal load and the dotted lines represent the cases without thermal 
load; in Fig. 25(b), each surface represents the intersection area variation with respect to Ny 
and Nx under certain thermal load. 



(a) 


(b) 


Figure 25: Variation of cross-section area with respect to Nx, Ny, and AT. 
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Actually, from those three basic numerical tests, the CODs can be expressed explicitly in terms 
of Nx, Ny, and A T. 


u = 0.592e _11 A/ x + 0.446e~ u N y - 1 ,082e 8 A T (22) 

v = 0.204 e u N x + '\A56e u N y - 1 ,036e 8 A7 (23) 

The corresponding interlaminar cross-sectional area AA is given by 

A=4xuxv 

= 0.121e 22 A/ 2 +0.516e 22 A/ 2 + 0.1 12e 15 AT 2 (24) 

+ 0.775e' 22 A/ x A/ y -0.834e“ 19 A/ x Ar-0.171e“ 18 A/ y AT 

From Eq. (24), it is obvious that the intersection area shows parabolic variations with respect to 
Nx, Ny, and AT, which can be seen clearly in Fig. 25. Equations (22) to (24), however, are only 

valid for the current model. If any condition of the model, such as delamination length, material 

properties and ply thickness, is changed, we have to redo the computations of those three basic 
cases. 

Effects of Delamination Shape and Delamination Length - The effects of delamination shape 
and the delamination length on the intersection area are investigated in this section. In previous 
section, the delamination fronts developed from the surface ply and the middle ply were 
connected by a pair of straight lines forming a rectangular delamination front. However, if three- 
dimensional J-integrals are computed along the crack front as shown in Fig. 26, it can be seen 
that the highest value occurs at the corner, i.e., the crack will first propagate at the corner. 
Hence, the delamination with a circular front might be a more appropriate scenario. 

Figure 26(a) shows the top view of quarter of the RVE with delaminations and Fig. 26(b) shows 
the distribution of J-integrals with different connections. The model is subjected to thermal load 
only. As we can see J-integral values remains as constant at the straight crack front and 
increases near the connecting corner. For the case of rectangular front, the J-integral value 
computed at corner, which is extremely high due to the discontinuity of first order derivative, 
may not be meaningful and is eliminated from Fig. 26(b). As we expected, the circular 
delamination front reduces the highest value of J-integral at corner and the bigger is the radius, 
the lower is the value. 
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Figure 26: J-integral distribution along crack front with different connect radius. 


The intersection areas in models with different connecting radius are computed. Figure 27 
shows the effects of the radius of the connecting circle. In all three basic cases, the cross- 
sectional area increases quadratically with respect to the radius. 



Figure 27: Effect of connecting radius of delamination on intersection area. 
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Moreover, the effect of delamination length is also investigated, as shown in Fig. 28. In this 
case, the rectangular front is still used and the delamination length changed from half of the ply 
thickness to twice the ply thickness. It can be seen that the intersection area is more sensitive to 
the delamination length than to the circle radius. It does not show the tendency to be a plateau 
as mentioned by Roy and Benjamin (2004). 



0 0.5 1 1.5 2 

Delamination length {a /h ) 


Figure 28: Effect of delamination length on intersection area. 

Single Delamination Model — Another phenomenon shown in Fig. 26(b) is that the steady state 
values on different sides are quite different even though both directions are under the same 
loading condition. Hence, delaminating from surface ply (high J-integral) is much easier to 
happen than delaminating from middle ply (low J-integral). From different view of the typical 
deformation of unit cell as shown in Fig. 29, it can be seen that delamination from middle ply 
(Fig. 29(a)) tries to close while the delamination from surface ply (Fig. 29(b)) tries to open up. 
Therefore, it is much possible that delamination only emanate from surface plies. Even if only 
one-side delamination appears, the path for gas leakage can still be formed as shown in Fig. 30. 
We call this model as “Single Delamination Model (SDM)”, and, correspondingly, the previous 
model is called “Double Delamination Model (DDM)”. 

In this case, the intersection area is smaller than the one in the previous model. Figure 31 
shows the comparison between the intersection areas computed using SDM and DDM under 
the same loading condition (N y /N x =0.5 and AT =-223 °C). The intersection area A computed 
using SDM also shows a parabolic variation, but the magnitude is only about half of that 
computed using DDM. 
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Figure 29: Typical deformation of unit cell from different views. 



Figure 30: Gas leakage path formed by the delamination from surface ply only. 



Figure 31 : Comparison between two different scenarios. 
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Effect of Lay-ups — Until now, only three-ply laminates are investigated. In reality, there are 
various lay-up possibilities. Choi (2005) did permeability tests for different laminates and found 
that for the same thickness the specimen with layers grouped together has higher permeability 
than the one with layers dispersed, i.e., laminate [0/90/0 2 /90/0] performs better than laminate 
[0 2 /90 2 /0 2 ]. In this section, the current modeling technique will be used to demonstrate this 
phenomenon. In this investigation, exact same lay-up and material properties are used as in 
Choi’s tests (Choi, 2005). For laminate [0 2 /90 2 /0 2 ], the unit cell is kept same as the one 
described in previous sections. For laminate [0/90/0 2 /90/0], in order to describe the unit cell, the 
micro-cracks in the 0°-plies and in 90°-plies are assumed to appear in same position and the 
crack densities in all 0°-plies or in all 90°-plies are equal, as shown in Fig. 32. Due to the 
symmetry of the model, only one-eighth of the unit cell is considered. Figure 32(b) shows the 
typical deformation of the unit cell. 

The crack densities in laminate [0 2 /90 2 /0 2 ] and in laminate [0/90/0 2 /90/0] are kept as same. The 
permeabilities in both laminates are computed using Eq. (11). Since the material constant C is 
to be determined, the permeability is normalized by the constant C. In Table 2 the normalized 
permeability k/C for both specimens are tabulated and compared with the experimental results 
by Choi (2005). From this table, it can be seen clearly that the grouped specimen has a much 
higher permeability than the dispersed specimen. Note that the ratios between the permeability 
of the two specimens are quite different in the models and the experiments. This is due to the 
fact that the modeling procedure does not precisely reflect the testing procedure. In Choi’s tests, 
the specimens are tested at room temperature after certain cryogenic cycles (cooling specimen 
down to cryogenic temperature and then warming it up to room temperature is defined as one 
cryogenic cycle) without mechanical loading. While in our modeling, crack density and 
delamination length is assumed and bi-axial loading is applied. Hence, the two sets of results 
can only be compared qualitatively for now. 
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Table 2: Comparison between Computed and Tested Permeability 


Laminates 

Normalized Permeability 
kIC 

Test results bv Choi (2005) 
(mol/sec/m/Pa) 

[0 2 /90 2 /0 2 ] 

1.26x10' s 

2.23x10 18 

[0/90/0 2 /90/0] 

3.30x1 O' 9 

1.16x10' 20 


Temperature Dependent Properties — In this section, the effect of the temperature dependent 
properties is investigated. In previous sections all material properties are assumed to be 
temperature independent. For graphite/epoxy composite laminates, the matrix properties largely 
depend on the temperature while the fiber properties do not. Hence, the composite properties 
that depend mainly on matrix properties, such as transverse modulus E 2 and shear modulus 
G 12 , vary remarkably with respect to temperature. Schulz (2005) and Speriatu (2005) have 
performed some experiments to measure the temperature dependent properties for 
graphite/epoxy composite IM7/977-2. Transverse modulus E 2 , shear modulus G 12 , and 
coefficient of thermal expansion a are measured at different temperature and fitted with 
polynomial functions. In the following simulation, these experimental results are used. 
Longitudinal modulus E-\ and Poisson’s ratios are assumed to be temperature independent, and 
shear modulus G 23 is assumed to follow the relation 


G 


23 


E 2 

2(l + o 2Z ) 


(25) 


In the simulation, constant strain boundary conditions are applied. Strains in both x-direction 
and y-direction are 1%. Thermal load is applied from room temperature to cryogenic 
temperature. Figure 33 shows the computed intersection area. The line with square symbols is 
the result with temperature independent material properties and the line with diamond symbols 
is the result with temperature dependent properties. From this figure, it can be seen that at 
cryogenic temperature the result with temperature dependent properties is almost 15% less 
than that with temperature independent properties. Hence, if temperature independent material 
properties were used, the simulation will give the upper bound of the permeability. 
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Figure 33: Effect of temperature dependent material properties on the intersection area. 

Progressive Permeability Prediction 

In this section, progressive gas permeability in laminates [0/90 n /0] (n = 1, 2, 4) is predicted in 
terms of applied load. Previously, the crack densities in laminates [0/90 n /0] are predicted under 
different biaxial loadings. A typical unit cell is taken to compute the intersection area. To keep 
same as the unit cell in crack density prediction, the surface cracks in the unit cell are not 
symmetric. The crack densities in this unit cell is taken as A x = A y = 2 cm' 1 . Delaminations from 
both surface ply and middle ply are introduced in this model. Due to symmetry, only one-fourth 
of unit cell is taken for analysis, as shown in Fig. 34. In this section, only mechanical loads are 
considered, so only two basic cases are needed for each unit cell: ( e x = 1%, e y = 0) and ( e x = 0, 
£ y = 1%). Figure 35 shows the deformation of the unit cell under loading case (ex = 1%, ey = 0). 
The deformation is enlarged by a factor of 10. 

The two basic cases of intersection areas for laminates [0/90„/0] are tabulated in Table 3. From 
this table, it can be seen that in both cases the intersection area increase as the thickness of the 
middle ply increase. A MATLAB program is then developed to predict progressive permeability 
and to plot the permeability versus the applied load. In this program, the predicted crack 
densities shown in Fig. 14 to Fig. 22 are fitted with cubic polynomials and then the intersection 
areas are computed from the two basic cases shown in Table 3. Finally, the progressive 
permeability is computed using Eq. (11). 

Figures 36 through 38 show the predicted permeability as a function of average stress a y . Each 
figure represents one loading ratio. In Fig. 36, N x IN y = 0.5, in Fig. 37, N x /N y = 1.0, and in Fig. 38, 
N x /N y = 2.0. From these figures, it can be seen that the initiation stresses of gas permeability 
increase as the thickness of middle ply decrease. After initiation, however, the gas permeability 
in the laminate with thinner middle ply increases faster than in the one with thicker middle ply. 
This follows the same trend as the evolution of crack density. Even though the intersection area 
in [0/90 4 /0] is the largest, the permeability in [0/90 4 /0] increases slowest comparing to the other 
laminates. Hence, it can be conclude that the permeability is dependent more on the crack 
density than on the intersection area. 
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For laminate [0/90/0], no permeability is predicted for N x IN y = 0.5 and N x IN y =1.0 under 
available stress level because the middle ply remains undamaged (see Fig. 14 and Fig. 17). 


Table 3: Typical Intersection Areas (m 2 ) for Laminates [0/90 n /0] 


Lay-ups 

[0/90/0] 

[0/90 2 /0] 

[0/90 4 /0] 

Case 1: £ x = 1 %, e y = 0 

8.32X10' 11 

1 0.5x1 0" 11 

12.5x1 0" 11 

Case II: e x = 0, e y = 
1% 

7.72X10' 11 

10.1X10' 11 

11.8X10' 11 
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Figure 35: Typical deformation of the unit cell of [0/90/0], 



Figure 36: Progressive permeability in laminates with loading ratio a = 0.5. 
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Figure 37: Progressive permeability in laminates with loading ratio a = 1 .0. 



Figure 38: Progressive permeability in laminates with loading ratio a = 2.0. 

Permeability Testing of Composite Laminates 

The permeability is defined by the amount of gas that passes through a given material of unit 
area and unit thickness under a unit pressure gradient in unit time. The SI unit of the 
permeability is mol/sec/m/Pa. Experiments were performed to investigate the effect of cryogenic 
cycling on permeability of laminated composites and to provide useful comparison of 
permeabilities of various composite material systems. 

Standard Test Method for Determining Gas Permeability 

The standard test method for determining gas permeability is documented in ASTM D1 434-82 
(Re-approved in 1997) “Standard Test Method for Determining Gas Permeability Characteristic 
of Plastic Film and Sheeting [ASTM].” The permeability can be measured by two methods, 
manometric determination method and volumetric determination methods. The permeability 
experiment using the monometric determination method is shown in Figure 39. The lower 
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pressure chamber is initially vacuumed and the transmission of the gas through the test 
specimen is indicated by an increase in pressure. 

The permeability is measured using volumetric determination as shown in Figure 40. The lower 
pressure chamber is exposed to atmospheric pressure and the transmission of the gas through 
the test specimen is indicated by a change in volume. The gas volume-flow rate is measured by 
recording the rise of liquid indicator in a capillary tube per unit time. The gas transmission rate 
(GTR) is calculated using the ideal gas law. The permeance is determined as the gas 
transmission rate per pressure differential across a specimen. The permeability is determined 
by multiplying permeance by the specimen thickness. 



The monometric determination method was not considered for this study since it is dangerous to 
handle the mercury compound, which is considered as a health hazardous material. Therefore, 
the permeability facility was constructed based on the volumetric determination method as 
shown in Figure 40. 
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Figure 40. Permeability experimental setup for volumetric determination method 


Permeability Apparatus 

This apparatus basically consists of a specimen placed between two chambers as shown in 
Figure 41. The test gas is pressurized in the upstream chamber. The gas permeates through 
one side of the specimen and escapes out of the other side. The escaped gas is collected in the 
downstream chamber and flows into a glass capillary tube. The amount of gas escaping per unit 
time is measured. The permeance is determined by gas transmitted rate and the partial 

pressure differential across the specimen. The permeability P is defined by the product of 
permeance P and the specimen thickness t. 

The gauge pressure of the gas in the upper upstream chamber is measured using a pressure 
transducer (P-303A from the Omega Engineering Inc.). The ambient pressure is measured by a 
barometric sensor (2113A from the Pasco Scientific). A precision pressure regulator provides 
constant gas pressure to the upstream chamber. The ambient temperature was measured using 
a glass capillary thermometer. 
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Figure 41 . Permeability testing apparatus. 


The specimens are mounted horizontally between the upstream and downstream chambers and 
clamped firmly by applying a compressive load as shown in Figure 42. The specimen is sealed 
with a rubber gasket and an O-Ring (38 mm inner diameter). A force gauge mounted at the top 
measures the compressive load to ensure that the same amount of compressive load is applied 
on the specimens for every test. The compressive load should be enough to prevent gas 
leakage, but should not damage the specimens. The upstream chamber has an inlet vent and 
an outlet vent. The inlet vent allows the gas flow into the upstream chamber and the outlet vents 
is used to purge the test gas to atmosphere (Figure 4-2). The downstream chamber has two 
outlet vents. One is used to purge the test gas to atmosphere and the other allows the gas flow 
to the glass capillary tube for measurements. The sensitivity of permeability measurement can 
be improved by increasing the gas transmitting area of a specimen. 
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The glass capillary tube is mounted on a rigid aluminum base horizontally to minimize the 
gravity effect on the capillary indicator and for easy reading of the scale marks on the capillary 
tube. Nettles (2001) found that there was no significant difference in the volumetric flow rate 
when capillary tube is placed vertically or slanted. The inner diameter of the glass capillary tube 
is 1.05 mm and the length is 100 mm. A magnifying glass is used to read the scale marks at the 
top of the meniscus of the liquid indicator. 

The liquid indicator in the glass capillary tube is used to measure the rate of rise of the liquid 
indicator. The rate is used to calculate the volume flow rate of the escaped gas across the 
specimen. Nettles (2001) investigated the effects on the volume flow rates by using various 
types of liquid. The volume flow rates obtained using water, alcohol, and alcohol with PhotoFlo® 
were not significantly different. In the present study, methyl alcohol is chosen as the liquid 
indicator since alcohol has low viscosity and weight. The methyl alcohol is colored with a blue 
dye to obtain precise readings on the scale marks. 

The primary investigation is the hydrogen permeability of a liquid hydrogen composite storage 
system. Since hydrogen gas is highly flammable and explosive when it mixed with air, it needed 
to be handled carefully during the test. Hence, other permeate gases were considered as a 
substitute for the hydrogen gas. The molecular diameter of various gases is listed in Table 4 
below (Nettles, 1995). To choose a permeate gas, the molecular diameter determined from 
viscosity measurement is mainly considered since the permeability is related with the volumetric 
flow rate directly. Since helium has the smallest molecular diameter, the helium predicts the 
permeation results higher than other gases. Therefore, in the study, helium was chosen as the 
permeate gas instead of hydrogen. 

Table 4: Molecular diameter of various gases 
(CRC Handbook of Chemistry and Physics, 54 th Edition) 


Type of Gas 

From Viscosity 

Molecular Diameter (cm) 

From van der 
Waal’s Equation 

From Heat 
Conductivity 

Helium 

1.9x1 O' 8 

2.6x1 O' 8 

2.3x1 0‘ 8 

Hydrogen 

2.4x1 O' 8 

2.3x1 O' 8 

2.3x1 O' 8 

Nitrogen 

3. 1x1 O' 8 

3.1x10' 8 

3.5x1 O' 8 


Specimen Description 

The permeability tests were performed with various composite material systems. The details of 
specimens are described in Table 5. The specimens Cl, C2, and C3 are various graphite/epoxy 
laminated composites. The specimen T1 is a textile composite. The specimen N1 is a laminated 
composite embedded with 36 pm aluminum oxide (Al 2 0 3 -alumina) nano-particles. The 
aluminum oxide was dissolved in alcohol and the compound was applied on a surface on a 
laminated prepreg using a paint blush. The nano-particles are capable of relieving the thermal 
stresses due to contraction of fiber and matrix phases. In addition, it can prevent the crack 
propagation in matrix phase by relieving the stress concentration at the crack-tip. The 
graphite/epoxy prepregs were fabricated with designed stacking orientations. The 
graphite/epoxy prepregs were cured in an autoclave and machined by a diamond saw to 
3x3 inch specimens at low speed to avoid fiber shattering. The surface is cleaned with acetone 
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and prepared carefully to avoid contamination or damages on the surface of the specimen 
during machining and subsequent handling. 

Table 5: Description of composite specimens 


Specimen 

Preprag Type 

Stacking Sequence 

Thickness (mm) 

Cl 

Graphite/Epoxy laminates 

[0/90/0/90/0/90] s 

1.52 

C2 

Graphite/Epoxy laminates 

[0 2 /90 2 /02]t 

0.787 

C3 

Graphite/Epoxy laminates 

[0/90/0 2 /90/0] t 

0.914 

C5 

Graphite/Epoxy laminates 

[0/90 2 /0] t 

0.533 

T1 

Plain weave textile(SP 
Systems SE-84) 

4 layers 

0.686 

N1 

Graphite/Epoxy laminates 
with nano-particles 

[0/90/NP/90/0] t 

0.483 


The specimens were subjected to cryogenic cycling at specified number of times, representing 
multiple refueling process of a space vehicle. A single cryogenic cycle consisted of cooling 
down from room temperature to cryogenic temperature and then warming up to room 
temperature. Initially, specimens were placed at room temperature (T = 300K). And, then the 
specimens completely submerged in an insulated container filled with liquid nitrogen. The 
specimens stayed in the container for approximately 2 minutes. When the specimen reached 
the boiling temperature of liquid nitrogen (T = 77K), the liquid nitrogen boiling disappears. The 
specimens were taken out of the container and placed at room temperature. The specimen was 
held at room temperature for approximately 5 minutes. The cryogenic cycling procedure is 
repeated for desired number of times. 

Testing Procedure 

Before starting the permeability test, a thin coat of silicon grease was applied on the gasket and 
an O-ring was placed to improve sealing of contact surfaces on a specimen. The excessive 
silicon grease was wiped out to avoid obstructions of permeated gas on the transmitting surface 
of the specimen. All outlet valves were opened initially to avoid sudden pressurization of test 
gas. The specimen was placed horizontally on the gasket of the upstream chamber. And then, 
the downstream chamber was placed on the top surface of a specimen. The specimen was 
mounted between the chambers. Both chambers were aligned and mated as close as possible. 
The specimen was mounted between the chambers and clamped firmly with a compressive 
force. Then, even distributed forces were applied to the sealing materials of the chambers. 

The test gas was admitted to the upstream chamber by opening the gas release valve of the 
gas tank. While all outlet valves remained opened, the test gas was filled in the upstream 
chamber and ventilated though the outlet vent to atmosphere. Any residual air in the upstream 
chamber was purged for 1 minute. The outlet valve on the upstream chamber was closed and 
the test gas was allowed to permeate across the specimen for a sufficient time to purge any 
residual air at downstream chamber. At this time, only test gas filled the chambers. When the 
outlet valve of the upstream chamber was closed, the upstream pressure increased slowly. The 
upstream pressure can be adjusted by controlling the pressure regulator. Sufficient time was 
allowed for attaining steady state of moving rate of the liquid indicator before beginning to take 
readings. The distance of rise of the liquid indicator was measured while the ambient pressures 
were recorded. 
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Calculations 

The volumetric methodology is used to calculate the permeability by measuring gas volume 
transmitted through a specimen. The rate of rise of the liquid indicator is measured using the 
capillary tube. The volume flow rate is calculated as follows. 

V r = slope x a c (26) 

where a c is the cross-sectional area of a capillary tube. The slope is the rate of rise of the liquid 
indicator in the capillary tube. The gas transmitted rate (GTR) is calculated using the ideal gas 
law as follows. 


GTR = 


PoX 

ART 


(27) 


where p 0 is ambient pressure, A is transmitting area of a specimen, R is the universal gas 
constant (8.3143x10 3 L-Pa/{mol K}) and T is ambient temperature. The permeance P is defined 
by the ratio between the gas transmission rate and pressure differential across the thickness of 
a specimen. Therefore, the permeance P is calculated as follows. 


P = 


GTR 

P-Po 


where p is the upstream pressure. The S.l. units of permeance are [mol/(m 2 -sPa)]. 


(28) 


According to the standard test method, the permeability P is defined by the product of 
permeance P and the specimen thickness t and “the permeability is meaningful only for 
homogeneous materials [ASTM].” In this study, although the laminated composites are 

considered as orthotropic materials, its permeability P is calculated by following the 
corresponding definition for isotropic materials [Nettles, 2001]. 

Calibration 

The position of capillary indicator is very sensitively to even small variations of testing conditions 
such as ambient pressure and temperatures caused by closing doors or air-conditioning system 
in the laboratory. For example, during one test, the variation of the ambient pressure was 
approximately 0.3% for a 13-hour period (Figure 43). For the permeability calculation, the 
ambient pressure is assumed to be constant. Flowever, the variation is large enough to cause 
error in measuring the rate of rise of the capillary indicator. The moving distance is needed to be 
calibrated to compensate for error due to the variation of ambient pressure. 
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Figure 43. Variation of ambient pressure for 13 hours at test condition 

The correction factor due to the variation of ambient pressure was calculated by performing the 
permeability test without applying the upstream pressure. An aluminum plate is used for the test 
to ensure that there is no gas permeation to the downstream chamber. Since there was not 
much variation in the ambient temperature, only ambient pressure causes changes in the 
position of the capillary indicator. After the outlet valve of the downstream chamber was closed, 
the displaced position of the capillary indicator and ambient pressure were measured 
simultaneously as shown in Figure 44. 



Figure 44. Variation of barometric pressure and indicator position as a function of time 
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The correction factor k is calculated by k = Ah /Ap where h is the moving distance of the 
capillary indicator and p is the ambient pressure. The average correction factor is found to be 
0.21 mm/Pa. Therefore, the corrected moving distance of the capillary indicator is calculated as 
follows: 


h corrected = h actual + k ■ Ap (29) 

where Ap is differential of ambient pressure between the beginning and end of the test. 

Permeability Test Results 

The permeability tests were performed with various composite material systems at room 
temperature. The permeability was measured at six different levels of upstream pressure. The 

average permeance P and average permeability P are tabulated in Table 6 for laminated 
composites, the results for textile composites are in Table 7 and the results for laminated 
composites with embedded with nano-particles in Table 8. 


Table 6. Permeability of laminated composites for various number of cryogenic cycles 


Specimen 

Cryogenic 

cycles 

Permeance, P 
(mol/sec/m 2 /Pa) 

Permeability, P 
(mol/sec/m/Pa) 

Logarithm 
of P 

Cl 

0 

5.60x10' 18 

8.54x1 O' 21 

- 20.1 


1 

1 ,52x10' 17 

2.32x1 O' 20 

-19.6 


5 

2.39x10' 17 

3.65x1 O' 20 

-19.4 


10 

2.39x10' 17 

3.65x1 O' 20 

-19.4 


20 

2.1 IxlO' 17 

3.22x1 O' 20 

-19.5 

C2 

0 

7.02x10' 18 

1.07x1 O' 20 

-20.0 


1 

1 ,06x10‘ 17 

1 .62x1 O' 20 

-19.8 


5 

1.47x1 0‘ 15 

2.23x10' 18 

-17.7 


10 

1 ,42x10" 15 

2.16x10' 18 

-17.7 


20 

1 ,49x10' 15 

2.27x10' 18 

-17.6 

C3 

0 

6.22x10' 18 

9.48x1 O' 21 

-20.0 


1 

7.56x10' 18 

1 ,15x10' 20 

-19.9 


5 

7.60x10' 18 

1 ,16x10' 20 

-19.9 


10 

8.37x10' 18 

1 .28x1 O' 20 

-19.9 

C5 

0 

5.85x10' 18 

8.92x1 O' 21 

-20.0 


1 

9.52x10' 16 

1 ,45x10' 18 

-17.8 


5 

8.67x10' 16 

00 

b 

X 

C\J 

00 

-17.9 


20 

CD 

b 

X 

00 

00 

1 ,34x10' 18 

-17.9 


NASA/CR— 2008-21 5440/PART3 


68 



Table 7. Permeability of textile composites for various number of cryogenic cycles 


Specimen 

Cryogenic 

cycles 

Permeance, P 
(mol/sec/m 2 /Pa) 

Permeability, P 
(mol/sec/m/Pa) 

Logarithm 
of P 

T1 

0 

4.79x1 O' 18 

7.30x1 O' 21 

-20.1 


1 

6.77x10" 18 

1.03x1 O' 20 

-20.0 


5 

8.41x10' 18 

1.28x1 O' 20 

-19.9 


20 

8.75x1 0‘ 18 

1.33x1 O' 20 

-19.9 


Table 8. Permeability of laminated composites embedded with nano-particles for various 

number of cryogenic cycles 


Specimen 

Cryogenic 

cycles 

Permeance, P 
(mol/sec/m 2 /Pa) 

Permeability, P 
(mol/sec/m/Pa) 

Logarithm 
of P 

N1 

0 

6.82x10' 18 

1.04x1 O' 20 

-20.0 


1 

2.72x10' 15 

4.15x10' 18 

-17.4 


5 

1 ,30x10‘ 14 

1 ,98x10' 17 

-16.7 


20 

9.83x10' 15 

1 ,50x10' 17 

-16.8 


The test results show the permeability increases as the number of cryogenic cycles increases 
(see Figure 45). The permeability increased rapidly and becomes constant with further increase 
of cryogenic cycles. For specimens C2 and C3, which have approximately same thickness, the 
permeability of the specimen C3 was lower since the specimen C3 has the plies dispersed and 
not grouped together compared to the specimen C2. 

The textile specimen T1 maintained constant permeability with the increase of cryogenic cycles. 
The textile composites resulted lower permeability than the laminated composites. The 
specimen N1 has same layer stacking orientations with the specimen C5 and nano-particles 
were dispersed between two 90-degree layers. Before cryogenic cycling, the permeabilities of 
the specimens N1 and C5 were approximately the same. However, as the number of cryogenic 
cycles increased, the permeability of the specimen N1 became higher. The use of nano- 
particles did not improve the permeability of laminated composites. 
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Figure 45. Logarithm of the permeability for composite specimens with increase of cryogenic 
cycles. 

Optical Microscopic Analysis 

The optical microscopic inspection was performed to evaluate the microcrack propagation and 
void content of various composite systems after cryogenic cycling. In the previous section, the 
experimental results showed that the permeability increases as the composite specimen 
underwent more cryogenic cycles. As the crack density increases, gas flow becomes easier 
though the specimen. Therefore, the microcrack propagation is correlated with cryogenic cycling. 

The specimens were cut through the center using a diamond saw. A LECO grinder/polisher was 
used for the sample preparation process. The rough edge through the center was ground with 
600-grit sand paper with water for 30 seconds. The fine grinding was performed with the 1000- 
grit and 1500-grit papers for 30 seconds. The surface of the edge was polished with the 58 pm 
aluminum oxide powder (Al 2 0 3 -alumina) dissolved in distilled water. The purpose of the 
lubricant is to both dissipate the heat from polishing and to act as a carrier for the abrasive 
materials. The ultrasonic cleaner was used to remove any abrasive particles and contaminants 
on the polished surface of a specimen. The optical analysis was conducted with a NIKON 
EPIPOT microscope. 

The laminated composite specimen C2 and the textile composite specimen T1 were chosen for 
optical inspection. The specimen details are described in Table 5. The microscopic images for 
the specimens were compared before and after cryogenic cycling. 
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(a) 



(b) 

Figure 46. Cross sectional view of the graphite/epoxy composite specimen C2 before cryogenic 
cycling: (a) 10X magnification; (b) 40X magnification. 

For the graphite/epoxy specimen C2 before cryogenic cycling, no microcrack propagation 
observed (see Figure 46). Some voids formed in the middle of the 90-degree layers and 
between the 0-degree layer and the 90-degree layers. The voids probably formed during 
composite fabrication in autoclave. When the graphite/epoxy prepreg was cured at high 
temperature, some air bubbles could have been trapped between layers. It is possible that 
external pressure applied on the vacuum bag was insufficient to remove the bubbles. 

After cryogenic cycling on the specimen C2, microcracks were observed in the 90-degree layer 
as shown in Figure 47. The fiber breaks were not observed in the in 0-degree layer since the 
thermal stresses were not large enough. The delaminations propagated along the middle of the 
90-degree layer where some voids were found. The transverse cracks branched with the 
delaminations. Since the gas can be transmitted through the transverse cracks across the 
specimen, the permeability increased after cryogenic cycling. 
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0.05 mm 


(b) 

Figure 47. Microcrack propagation on the graphite/epoxy composite specimen C2 after 
cryogenic cycling: (a) 10X magnification; (b) 40X magnification. 

For the textile composite specimen T1 before cryogenic cycling, no microcracks were observed 
as shown in Figure 48. Voids were observed at the location where two yarns are merging. 



Figure 48. Cross sectional view of the textile specimen T1 before cryogenic cycling, 10X 
magnification. 
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(b) 

Figure 49. Microcrack propagation in the textile specimen 
T1 after cryogenic cycling: (a) 10X magnification; 

(b) 40X magnification. 

After the textile specimen T1 underwent cryogenic cycles, microcrack were observed in 90- 
degree yarn as shown in Figure 49. The microcracks propagated in the in-plane direction of the 
specimen. Since transverse cracks were not propagated, the permeability of textile composites 
was almost the same before and after cryogenic cycling. The transverse cracks and 
delaminations provide the leakage path through composite laminates and thus directly related to 
the permeability. For the laminated composite specimen, the microscopic results showed that 
transverse crack propagation and delaminations of composite laminates. For the textile 
composite specimen, the transverse cracks propagation was not observed, which resulted low 
permeability of textile composites even after cryogenic cycling. 

Conclusions 

In this work, a permeation model for composite laminate is developed based on Darcy’s Law. 
The permeability constant is a function of crack density, intersection area and a material 
constant. From 2-D FE analysis, the problem is decoupled into two parts: crack density 
prediction and intersection area computation. The crack density in cross-ply laminate under 
biaxial loading is predicted as a function of applied load using an energy-release-rate-method 
based on 3-D FE analysis. Then the intersection area for a typical unit cell is computed using 
3-D FEM because of the significant local effect. The effects of various factors on the intersection 
area are investigated, and it has been found that the intersection area is not only dependent on 
applied load but also dependent on delamination length, shape of the delamination front, 
temperature dependent material property and lay-ups. The progressive permeability in cross-ply 


NASA/CR— 2008-21 5440/PART3 


73 




laminate with different middle ply is predicted and it has been found that the initiation stress of 
leakage and the increase rate of leakage are dependent on the thickness of the middle ply. The 
laminate with thinner middle ply has higher initiation stress, but also has higher increase rate 
once it start to leak. 

The experimental analysis was performed to measure the gas permeability of various composite 
material systems for the liquid hydrogen composite tank and the effect of cryogenic cycling of 
composite laminates on the permeability was investigated. The permeability test facility was 
constructed following the standard test method documented in ASTM D1 434-82. The 
permeability test was conducted at room temperature. The sensitivity of permeability 
measurement can be improved by increasing the gas transmission area of a specimen. Since 
hydrogen is flammable and explosive when it mixed with air, the helium gas is substitute for 
hydrogen. The calibration is performed to compensate the ambient pressure differences during 
the test. The correction factor is found as 0.21 mm/Pa. The actual moving distance of the liquid 
indicator is calibrated by the correction factor. 

The permeability of laminated composite was found to degrade after undergoing cryogenic 
cycling process. In thin laminate specimens the permeability increases significantly after 
cryogenic cycling. For thick laminated specimens the increase of permeability is less. The textile 
composite specimen has lower permeability than laminated specimens and the variation of 
permeability is very small with the increase of cryogenic cycles. The laminated composites were 
embedded with nano-particles, which are capable of reducing thermal stresses in the matrix 
phase. However, as results, the nano-particles did not show the improvement on permeability. 
Further studies are needed to investigate the effects of nano-particles. The optical analysis was 
performed to investigate the microcrack propagation and void contents of test specimens. The 
transverse cracks and delaminations provide the leakage path through composite laminates and 
thus directly related to the permeability. The optical micrographic analysis was performed to 
investigate the microcrack propagation and void contents of test specimens. For the laminated 
composite specimen, the microscopic results showed that transverse crack propagation and 
delaminations of composite laminates. For the textile composite specimen, the transverse 
cracks propagation was not observed, which resulted in low permeability of textile composites 
even after cryogenic cycling. 
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Abstract 

This work focuses on measuring and predicting the heat transfer rates associated with the 
complex, unsteady, nucleate flow boiling heat transfer that occurs during cryogenic chilldown. 
An inverse numerical technique is used to extract the heat transfer coefficient from time 
dependent measurements of wall temperature. This technique utilizes knowledge of the local 
temperature of the pipeline coupled with numerical simulations of the unsteady heat transfer 
through the wall of the pipeline to determine the heat transfer coefficients. The spatially varying 
heat transfer coefficients are averaged and compared against four nucleate flow boiling heat 
transfer correlations. None of the correlations are satisfactory in accounting for the spatial 
variations in heat transfer rate. 


Introduction 

The recent interest in space exploration has placed a renewed focus on rocket propulsion 
technology. Cryogenic propellants are the preferred fuel for rocket propulsion since they are 
more energetic and environmentally friendly other current storable fuels. Voracious boiling 
occurs while transferring these fluids through a pipeline that is initially in thermal equilibrium with 
environment. This phenomenon is referred to as line chilldown. Large temperature differences, 
rapid transients, pressure fluctuations and the transition from the film boiling to the nucleate 
boiling regime characterize the chilldown process. 

Although the existence of the chilldown phenomenon has been known for decades, the process 
is not well understood. Attempts have been made to model the chilldown process; however the 
results have been fair at best. A major shortcoming of these models is the use of heat transfer 
correlations that were developed for steady, non-cryogenic flows. The development of reliable 
heat transfer correlations for cryogenic chilldown has been hindered by the lack of experimental 
data. 

An experimental facility was constructed that allows the flow structure, the temperature history 
and the pressure history to be recorded during the line chilldown process. Details of the 
experimental facility may be found in the work of Velat [1] and Jackson et al. [2], The 
temperature history is then utilized in conjunction with an inverse heat conduction procedure 
that allows the unsteady heat transfer coefficient on the interior of the pipe wall to be extracted. 
It is observed that two predominant heat transfer regimes are encountered during chilldown. 
The first is the film boiling regime in which the liquid film is supported on a thin vapor layer 
adjacent to the wall, thus preventing the liquid from coming into contact with the interior pipe 
wall. This regime is immediately followed by the nucleate boiling regime in which the interior 
pipe wall temperature falls below the Leidenfrost point, hence the liquid contacts the inner wall 
of the pipe. This transition is marked by a sudden and significant decrease in wall temperature 
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that is directly related to an increase in heat transfer associated with nucleate boiling. The 
prediction of heat transfer rates in the film boiling regime has been considered in another work 
[2], This work explores cryogenic heat transfer in the nucleate flow boiling regime. The extracted 
heat transfer coefficients for the nucleate boiling regime are used to evaluate the predictive 
capability of various heat transfer correlations. 


Experimental Facility 

The experimental facility constructed to investigate the chilldown process is schematically 
illustrated in Fig. 1. The liquid nitrogen supply is stored in a 1580 kPa dewar. The driving 
potential for the nitrogen flow through the facility is supplied by the dewar pressure. Upon 
leaving the dewar the nitrogen is passed through a heat exchanger where it is cooled prior to 
entering the 12.5 mm I.D. vacuum jacketed Pyrex visual test section. The flow structure is 
captured using a high resolution CCD camera system. 



Figure 1. Cryogenic two-phase flow experimental facility. 


The flow then enters the heat transfer measurement section of the facility, located 70 cm 
downstream of the visual test section. In this section two inline thermocouples measure the 
temperature of the fluid at the inlet and outlet, while a series of thermocouples are placed 
circumferentially on the surface of the pipe and on the outer surface of the melamine insulation, 
so that the temperature profile on the top, bottom and either side of the pipe is obtained, see 
Fig 2. A throttling valve is place after the heat transfer section to allow the mass flow to be 
varied. The two-phase flow is then passed through two heaters, which vaporize any remaining 
liquid, to ensure single-phase vapor enters the venturi flow meter. A thermocouple and 
pressure transducer located prior to the venturi allow for the correction of compressibility effects 
to be taken into account in the mass flow rate computations. The vapor is then collected in an 
expansion tank, and evacuated from the facility via a ventilation system. 



Figure 2. Thermocouple arrangements along the pipe wall and along the insulation of the heat 
transfer section. 
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A Validyne differential pressure transducer measures the pressure drop across the venturi. The 
analog signals from the instruments are transferred and recorded by a 16-bit Measurement 
Computing data acquisition board and multiplexer system. The instruments have been 
independently calibrated and tested to ensure they function properly. The reported accuracy of 
the Validyne pressure transducers is ±0.25% of full scale (0-1380 kPa for the static pressure 
transducers, 0-8.62 kPa for the differential pressure transducer across the test section, and 
0-86.2 kPa for the differential pressure transducer across the venturi) that was confirmed via 
calibration. The static pressure transducers are accurate to ±3.5 kPa while the differential 
pressure transducer across the visual test section is accurate to within ±0.06 kPa and that 
measuring the pressure drop across the venturi is accurate to within ±0.35 kPa. Over the large 
temperature range experienced in these experiments the thermocouples used are accurate to 
within ±1.7 °C. 


Sample Measurements 

Experiments are carried out for mass fluxes that range from 66 to 500 kg/m 2 s; and vapor 
qualities vary from 0.004 to 1 . A typical plot of the temperature profile for the outer surface of the 
pipe in shown in Fig. 3; only one side wall temperature profile is shown since the difference 
between the to side thermocouples was within the error range of the thermocouples. The 
transition between the film boiling regime and the nucleate boiling regime is marked by the 
sudden and significant decrease in the wall temperature as highlighted in Fig. 3. The 
corresponding mass flux, vapor volume fraction and quality are also shown in Figs. 4, 5, and 6 
respectively. The mass flux and vapor volume fraction profiles highlight the inherent unsteady 
nature of the chilldown process. The vapor volume fraction was computed from the digital 
images of the flow structure and smoothed by finding the best smooth curve fit to the 
measurements; these values from the curve fit are then used to compute the vapor quality using 
the Zuber and Findlay correlation [3]; the result is illustrated in Fig. 6. 

This work focuses on the nucleate flow boiling regime and as such there are certain difficulties 
that are encountered while extracting the heat transfer coefficients from this region. The major 
difficulty results from the fact that the actual wetted perimeter of inner wall of the tube is not 
known and cannot be measured. In the film boiling regime the liquid film is supported on a thin 
layer of vapor adjacent to the pipe wall thus it is not in contact with the wall. Once the flow 
transitions to the nucleate boiling regime the liquid comes into contact with the wall, but this will 
not happen instantaneously over the entire region where liquid film is seen to reside. There will 
be a time lapse between the time the liquid at the bottom of the film touches the wall and the 
time the liquid at the top of the film touches the wall. Hence it is possible to overestimate the 
size of the wetted region when assuming that the wetted perimeter is equal to the liquid film 
height during the onset of nucleate flow boiling. The consequences of overestimating the size of 
the wetted region size are discussed when considering the inner surface boundary condition in 
the following section. 
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Figure 3. Outer wall temperature profile during chilldown process. Thermocouples are positioned 

at the top, bottom, left and right sides of the pipe. 



Figure 4. Mass flux variation with time for the chilldown process. 



Figure 5. Vapor volume fraction variation with time for the chilldown process. 



Figure 6. Vapor quality variation with time for the chilldown process. 
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Extracting the Heat Transfer Coefficient 

In order to extract the heat transfer coefficient from the temporal profile of wall temperature; an 
inverse procedure is employed which varies from that of traditional inverse heat conduction 
methods. It does not require a system of least-square equations to be solved. Some of these 
traditional methods are reviewed by Ozisik [4], 

The process used here for extracting the transient heat transfer coefficient on the inside of the 
pipe involves a number of iterative procedures, which are reported here as three major steps. 
These steps are illustrated in Fig. 7 and are carried out for each instance of time. These steps 
include (1 ) guess the heat transfer coefficient on the inside of the pipe, (2) knowing the heat 
transfer coefficient on the outside of the pipe via calibration, the temperature field is calculated, 
and (3) the temperatures calculated at the outer wall are then compared with those measured. If 
the temperatures match, the guessed heat transfer coefficient is taken as the actual heat 
transfer coefficient, otherwise a new guess is made and the process is repeated until the 
computed and measured temperatures match. 

Computing the Temperature Field in the Pipe Wall 

The unsteady 3-D form of the heat conduction equation in cylindrical coordinates, see Fig. 8, is 
written as follows: 
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Figure 7. Flow chart for transient heat transfer coefficient extraction. 

The temperature variations experienced by the pipe during the chilldown process are significant, 
and thus the thermal properties (k and c ) of the pipe material vary significantly. These 
variations are taken into account by assigning the thermal properties as a function of the 
temperature at any given point. 

A finite volume formulation is used to disceritize Eq. (2); a backward Euler scheme is employed 
for the temporal term and a central difference scheme is used for the spatial terms. The system 
of equations that result are solved using the Alternating Direction Implicit (ADI) method 
described in [5]. 



Figure 8. Coordinate system for heat conduction through the pipe wall. 


In order to solve the system of equations, the energy entering the system from the ambient must 
be known; hence the heat transfer coefficient for the outside surface of the pipe insulation must 
be determined. This is accomplished through a steady-state calibration process in which cool 
nitrogen vapor is passed through the heat transfer section at various mass flow rates. A 
relationship between heat flux through the pipe insulation versus temperature difference 
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(ambient temperature less the insulation surface temperature) is determined. A constant heat 
transfer coefficient is approximated as shown in Fig. 9. The slope of the line is the heat transfer 

coefficient for the outer pipe surface. The measurement of the heat flux into the pipe, q" , is 
quite difficult to measure due to the small temperature rise in the cryogenic vapor; thus there is 
some scatter in Fig. 9. Flowever, the measured heat transfer coefficient, 4.38 W/m 2 -K, is quite 
small, compared with convective heat transfer on the inner pipe wall. The errors in the 
calibration have a negligible impact on the computed thermal field in the pipe wall. 



Figure 9. Calibration for determining the outer pipe surface heat transfer coefficient. 


The boundary condition for the outer pipe surface is written, 
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which is non-dimensionalized to give, 
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Here h is the heat transfer coefficient, and subscripts out and » denote the outer surface of the 
pipe, and the outer insulation surface, respectively. The inner pipe boundary condition is written 
as, 
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Here the subscripts ; and fl denote the inner surface of the pipe and the fluid flowing through 
the pipe, respectively. The fluid temperature is determined using an internal thermocouple that 
measures the temperature at the inlet of the heat transfer section. 

The heat transfer coefficient on the inside of the pipe, h t , is the quantity that is guessed during 
the calculation. The two-phase flow structure present for much of the experiment is stratified, 
and there is a significant difference between the temperature at the top of the pipe and the 
temperature at the bottom the pipe. This must be due to circumferential variations in the heat 
transfer coefficient. This problem is dealt with adequately by dividing the interior surface of the 
pipe into three distinct sections within which the heat transfer coefficient is assumed constant. 
Region 1 is such that ()<«<«,, and k = h lop , region 2 is such that «,<«<«,, and h t =h siJe , and 

region 3 is such that a 2 <a<a 3 , and h i =h bottom (see Fig. 10). In the film boiling regime it is 
adequate to allow all three regions to be equal in size, hence, a x =a 2 -a x =n-a 2 . However in 
the nucleate boiling regime, where the temperature at the bottom changes suddenly by a 
significant amount very high heat transfer coefficients are present. Hence it is important not to 
overestimate the size of the region in which nucleate boiling occurs. If the region is 
overestimated then the amount of energy removed from the pipe would be too large in the 
adjacent region (in this case region 2) making it impossible to match the outer wall temperature. 
This is handled adequately by reducing the size of region 3 by increasing the value of a 2 . When 
necessary, it is sufficient to reduce the size of region 3 by a factor of 2. For modeling purposes, 
a more systematic approach to describing different regions during nucleate boiling is desirable. 
However, the present approach is sufficient to match the measured and computed wall 
temperature variations. 


i a 



Region 2 

hphside 


Figure 10. Assumed variation of heat transfer coefficient on the inside surface of the pipe. 
Iteration Process for Guessing the Inner Heat Transfer Coefficient 

Any number of methods may be used to obtain a new guess for the inner heat transfer 
coefficient; a systematic method for iterating is recommended. In this work after the initial guess 
is made subsequent guesses are determined using linear interpolation or linear extrapolation. 
Hence the new guess for inner heat transfer coefficient is given by, 
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where s is a scaling factor that reduces oscillations and takes the value 0.3 in this work. The 
subscripts pres , prev , new , and known denote the present, previous, new, and known quantities 
respectively. 

Test for Convergence 

The final step in the process is checking the computed temperature against the measured 
temperature. This is done using the temperatures at the top, side and bottom of the pipe. Once 
the computed temperature is within ±1x1 O' 8 K, the temperatures are considered a match. Care 
must be taken when determining the limit within which to consider the temperatures to be 
matched; if the value selected is too small the computation time increases significantly and the 
marginal improvement in accuracy does not justify the increased computational cost. 


Testing the Inverse Procedure 

In order to assess the performance and validity of this inverse method, it is first used to 
calculate the heat transfer coefficients for a single-phase flow simulation in which the heat 
transfer coefficient is known and varies in time. Second it is used to calculate the actual heat 
transfer coefficient for single-phase nitrogen gas flowing through the experimental facility, with 
the results being compared with the predictions of the Dittus-Boelter correlation for cooling. In 
the first test case, the heat transfer coefficient follows a parabolic path with time. A comparison 
of the specified heat transfer coefficient with that extracted using the inverse method is shown in 
Fig. 11. In Fig. 11 it is seen that a parabolic varying heat transfer coefficient is captured quite 
reliably with the inverse approach. 

A comparison of the extracted single-phase heat transfer coefficient with the Dittus-Beoelter 
correlation for flowing nitrogen gas is shown in Fig. 12. The comparison is quite good. The heat 
transfer coefficient is not constant with time since the mass flux fluctuates around a mean of 
196 kg/m 2 -s. 


The inverse technique is next applied to the experimental data for the chilldown process in the 
nucleate boiling regime. The unsteady heat transfer coefficients on the inner surface of the pipe 
are extracted over regions 1, 2, and 3. An average two-phase heat transfer coefficients is 
defined and compared with the correlations of Gungor and Winterton [7], Kandlikar [8], Muller- 
Steinhagen and Jamialahmadi [9] and Wojtan et al. [6,10], 
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Figure 11. Computation of a parabolic varying heat transfer coefficient using the inverse method. 
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Figure 12. Comparison of heat transfer coefficient computed using the inverse procedure and the 
Dittus-Boelter correlation for single-phase nitrogen gas flow. 

Computing the Two-Phase Heat Transfer Coefficient 

Since the size of the regions may change from one data set to the other, it is convenient to 
define an average two-phase heart transfer coefficient. This is done by integrating the heat 
transfer coefficients assigned on the inner surface over the perimeter of the inner surface. This 
approach is also employed by Wojtan et al. [6], Hence the average two-phase heat transfer 
coefficient, h lp , is computed as 


u -^L h , 2 ( g 2~«i) , 2^-^) 

n tp ~ 0 n top ^ 0 n side ^ 0 n l 

2 n 2 n 2 n 


( 8 ) 


This allows the extracted heat transfer coefficients to be compared with the existing correlations. 
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Results and Discussion 


Figure 13 illustrates the variation of the average two-phase heat transfer coefficient extracted 
from the experimental data with time for the chilldown process. It is clearly seen that the 
transition from the film boiling regime to the nucleate flow boiling regime is accompanied by an 
order of magnitude increase in the average two-phase heat transfer coefficient. The robust 
nature of the computational method is demonstrated in that it is able to handle the step change 
in heat transfer coefficient from the film to nucleate boiling regime. 



Figure 13. Average two-phase heat transfer coefficient variation with time. 

The nucleate flow boiling heat transfer coefficients extracted from the experimental data are 
presented in Table 2. The mass flux, vapor quality, flow regime, and saturation temperature are 
shown. In addition, the fraction of the perimeter occupied by each region, the average wall 
temperature, and the contribution to the two-phase heat transfer in each region are tabulated. It 
is observed that the influence of region 1 on the average two-phase heat transfer coefficient in 
the nucleate flow boiling regime is small when compared to the influence of regions 2 and 3. 
This is as a result of the stratified nature for the flow structure, in which the liquid resides in the 
lower regions of the pipe while the vapor resides in the upper region. By examining the 
contributions from each region to the total heat transfer, it is observed that in some instances 
region 2 has a larger contribution than region 3. This occurs because of the unsteady nature of 
the flow which results in instances where region 2 is periodically wetted resulting in a thinner 
liquid film than in region 3. Hence the heat transfer rate in region 2 will be greater as a result of 
the lower thermal resistance of the liquid film. 

In examining the fraction of the area assigned to each region it is evident that the flow structure 
influences the size of each region. For annular flow and stratified-wavy flow in which the liquid 
film height is greater than the radius of the pipe, the three regions may be constructed such that 
they are of equal size. While in the intermittent flow regime the size of each region is adjusted 
slightly. However, when the flow is in the stratified-wavy regime, a significant difference in the 
size of regions is observed. This results from the inability to determine the actual wetted area of 
the pipe. If the size of region 3 is overestimated the energy removed from the pipe wall would be 
too large in the adjacent region (in this case region 2) making it impossible to match the outer 
wall temperature. This problem is encountered mainly in the stratified-wavy regime, in which the 
height of the liquid film is less than the pipe radius. It is not encountered with the other flow 
structures since the wetted area is not overestimated. Provided the wetted region is not 
overestimated, the size of regions 1, 2, and 3 may have different values, and yet suitable 
agreement with measured data is found. Therefore, the solution presented is not unique. 
Nevertheless, it is found that the average heat transfer coefficient hardly changes with variations 
in the size of regions 1 , 2, and 3. In order to maintain consistency, these regions were taken to 
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be of equal size, unless the solution required differently. For this work the regions sizes were 
chosen to maintain as much symmetry as possible while ensuring that the size of the wetted 
region is not overestimated. The database of cryogenic nucleate flow boiling heat transfer 
coefficient for chilldown presented in Table 2 is the only one we are aware of. 

Comparisons of the computed average flow boiling heat transfer coefficients with those 
determined experimentally are shown in Figs. 14, 15, 16, and 17. Figs. 14, 15, 16, and 17 show 
comparisons with the Gungor and Winterton, Kandlikar, Miiller-Steinhagen and Jamialahmadi, 
and Wojtan et al. correlations, respectively. It is observed that the Miiller-Steinhagen and 
Jamialahmadi correlation performs the best, predicting 47.6% of the data within the 25% error 
band. This is followed by the Gungor and Winterton correlation with 38.1%, then the Wojtan et 
al correlation with 21.4% and finally the Kandlikar carrleation with 9.5%. 

The Kandlikar correlation and Gungor and Winterton correlation both over predict the data. Thus 
in applying these correlations care must be taken in order to account for the region that is 
actually wetted during the flow boiling process. The Wojtan et al. correlation under predicts the 
data. Although this correlation is the only one tested that takes into account the flow structure, it 
has not been extensively tested within the low quality range encountered in the chilldown 
experiments. The Miiller-Steinhagen and Jamialahmadi performs the best for the present data, 
however this correlation does not account for the flow structure. 

Further work is needed in developing a reliable correlation for nucleate flow boiling during 
cryogenic chilldown. 



Figure 14. Comparison of predicted and measured average nucleate flow boiling heat transfer 
coefficients using Gungor and Winterton correlation. 



Measured Heat T ransfer Coefficient (W/m 2 K) 

Figure 15. Comparison of predicted and measured average nucleate flow boiling heat transfer 

coefficients using Kandlikar correlation. 
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Figure 16. Comparison of predicted and measured average nucleate flow boiling heat transfer 
coefficients using Miiller-Steinhagen correlation. 



Figure 17. Comparison of predicted and measured average nucleate flow boiling heat transfer 

coefficients using Wojtan et al. correlation. 
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Table 2. Summary of measured nucleate flow boiling heat transfer coefficients using inverse method for 
regions 1 , 2, and 3. SW (0.5+) denotes stratified wavy flow where the height of the liquid film greater than 
djl ; I denotes intermittent flow; SW (0.5-) denotes stratified wavy flow where the height of the liquid film 

less than dj 2 ; A denotes annular flow. 


Flow 

Regime 

Mass 

Flux 

Vapor 

Quality 

T S at 

Fraction of 
area & Avg. 
Wall Temp for 
Region 1 

Fraction of area 
& Avg. Wall 
Temp for 
Region 2 

Fraction of 
area & Avg. 
Wall Temp for 
Region 3 

Heat Transfer 
Coefficient 
Contribution 
of Region 1 

Heat Transfer 
Coefficient 
Contribution 
of Region 2 

Heat Transfer 
Coefficient 
Contribution 
of Region 3 

ht P 


kg/m 2 s 


°K 

cx-i/ti 

°K 

(a2-ai)/7r 

°K 

[n-a2)ln 

°K 

W/m 2 K 

W/m 2 K 

W/m 2 K 

W/m 2 K 

SW (0.5+) 

144 

0.063 

98.7 

0.333 

151.5 

0.333 

131.0 

0.333 

105.3 

40 

336 

1828 

2203 

SW (0.5+) 

77 

0.016 

99.9 

0.333 

153.8 

0.333 

130.8 

0.333 

109.4 

25 

145 

1390 

1561 

SW (0.5+) 

78 

0.016 

99.9 

0.333 

151.2 

0.333 

125.2 

0.333 

110.1 

17 

382 

1454 

1852 

SW (0.5+) 

66 

0.006 

98.7 

0.333 

156.8 

0.333 

141.8 

0.333 

116.0 

18 

239 

1008 

1265 

SW (0.5+) 

66 

0.005 

99.2 

0.333 

155.2 

0.333 

136.6 

0.333 

106.9 

22 

458 

1731 

2211 

SW (0.5+) 

75 

0.004 

98.6 

0.333 

152.9 

0.333 

131.0 

0.333 

107.8 

24 

823 

1660 

2508 

SW (0.5+) 

234 

0.040 

101.8 

0.333 

151.0 

0.333 

127.3 

0.333 

112.1 

53 

654 

1174 

1881 

SW (0.5+) 

200 

0.042 

102.5 

0.333 

145.3 

0.333 

111.5 

0.333 

108.8 

56 

1543 

1289 

2888 

SW (0.5+) 

199 

0.034 

99.1 

0.333 

139.6 

0.333 

111.2 

0.333 

108.2 

146 

1503 

1138 

2787 

SW (0.5+) 

201 

0.036 

101.2 

0.333 

134.2 

0.333 

110.0 

0.333 

107.1 

180 

1224 

1025 

2429 

SW (0.5+) 

320 

0.011 

94.5 

0.333 

138.9 

0.333 

122.6 

0.333 

107.8 

51 

416 

1047 

1514 

SW (0.5+) 

310 

0.019 

94.4 

0.333 

135.3 

0.333 

107.7 

0.333 

102.2 

10 

1154 

1209 

2374 

SW (0.5+) 

302 

0.014 

94.2 

0.333 

132.6 

0.333 

105.3 

0.333 

103.2 

16 

1188 

384 

1588 

SW (0.5+) 

263 

0.014 

99.2 

0.333 

145.6 

0.333 

122.9 

0.333 

108.1 

34 

681 

1252 

1967 

SW (0.5+) 

262 

0.009 

99.7 

0.333 

142.0 

0.333 

111.4 

0.333 

107.5 

5 

1310 

1066 

2381 

1 

207 

0.008 

91.6 

0.417 

126.7 

0.292 

99.4 

0.292 

100.3 

2 

773 

1015 

1790 

1 

156 

0.008 

101.7 

0.417 

152.1 

0.292 

131.6 

0.292 

115.1 

32 

353 

915 

1300 

1 

203 

0.009 

102.2 

0.417 

148.6 

0.292 

118.1 

0.292 

111.3 

19 

986 

875 

1880 

1 

187 

0.009 

102.0 

0.417 

145.7 

0.292 

113.8 

0.292 

111.6 

19 

1199 

336 

1554 

1 

152 

0.009 

100.3 

0.417 

142.9 

0.292 

113.1 

0.292 

110.6 

33 

952 

359 

1345 

1 

182 

0.007 

98.7 

0.417 

140.1 

0.292 

112.3 

0.292 

108.4 

54 

820 

757 

1631 

1 

138 

0.007 

99.5 

0.417 

137.2 

0.292 

110.1 

0.292 

108.0 

49 

1322 

560 

1932 

SW (0.5-) 

94 

0.066 

83.8 

0.417 

121.0 

0.417 

107.2 

0.167 

93.1 

33 

311 

903 

1247 

SW (0.5-) 

128 

0.023 

85.9 

0.417 

119.5 

0.417 

93.0 

0.167 

90.2 

18 

1524 

3555 

5097 

SW (0.5-) 

90 

0.095 

86.7 

0.417 

117.9 

0.417 

101.0 

0.167 

89.2 

39 

645 

1165 

1848 

SW (0.5-) 

131 

0.089 

86.3 

0.417 

124.3 

0.417 

110.6 

0.167 

97.5 

37 

289 

777 

1103 

SW (0.5-) 

145 

0.091 

85.9 

0.417 

118.6 

0.417 

91.8 

0.167 

92.8 

21 

1215 

1019 

2255 

SW (0.5-) 

123 

0.038 

86.6 

0.417 

120.2 

0.417 

100.9 

0.167 

94.6 

19 

931 

1153 

2103 

SW (0.5-) 

135 

0.052 

84.7 

0.417 

118.1 

0.417 

105.4 

0.167 

97.2 

70 

523 

738 

1330 

SW (0.5-) 

141 

0.095 

84.7 

0.417 

115.4 

0.417 

90.4 

0.167 

94.4 

6 

2030 

1770 

3806 

SW (0.5-) 

181 

0.020 

85.4 

0.417 

115.5 

0.417 

99.4 

0.167 

96.3 

181 

993 

729 

1904 

SW (0.5-) 

173 

0.068 

87.1 

0.417 

119.7 

0.417 

96.4 

0.167 

96.0 

57 

1431 

775 

2264 

SW (0.5-) 

175 

0.069 

86.7 

0.417 

119.1 

0.417 

97.2 

0.167 

96.0 

45 

1310 

864 

2218 

SW (0.5-) 

89 

0.035 

91.8 

0.417 

140.1 

0.417 

120.5 

0.167 

101.3 

37 

402 

903 

1342 

SW (0.5-) 

116 

0.020 

93.8 

0.417 

138.9 

0.417 

114.7 

0.167 

103.7 

27 

376 

677 

1081 

SW (0.5-) 

254 

0.049 

91.8 

0.417 

124.0 

0.417 

107.1 

0.167 

105.5 

113 

1006 

586 

1706 

A 

345 

0.044 

94.5 

0.333 

130.8 

0.333 

114.7 

0.333 

108.4 

41 

676 

936 

1654 

A 

344 

0.044 

94.4 

0.333 

128.3 

0.333 

104.8 

0.333 

104.4 

3 

1304 

841 

2149 

A 

397 

0.081 

97.1 

0.333 

127.3 

0.333 

110.4 

0.333 

104.6 

107 

929 

1189 

2225 

A 

502 

0.129 

96.6 

0.333 

124.9 

0.333 

104.7 

0.333 

106.0 

42 

1664 

497 

2203 

A 

352 

0.039 

95.2 

0.333 

127.4 

0.333 

111.2 

0.333 

106.5 

49 

756 

877 

1682 

A 

415 

0.053 

95.0 

0.333 

124.9 

0.333 

104.0 

0.333 

105.7 

21 

1403 

383 

1808 
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Conclusion 


1 ) The lower portion of the pipe (regions 2 and 3) contributes the majority of the heat transfer in 
the nucleate flow boiling regime. 

2) The unsteady nature of the chilldown process may lead to region 2 being periodically wetted 
with a very thin liquid film and large heat transfer coefficient. 

3) The inverse technique has proved to be computationally robust in determining nucleate flow 
boiling heat transfer coefficients for the chilldown process. 

4) None of the nucleate flow boiling correlations tested is satisfactory, yet of those tested, the 
Miiller-Steinhagen and Jamialahmadi correlation has the best agreement with the data. 

Patents 

There was no patent application. 


Publications 

Velat, C., Jackson, J., Klausner, J.F., and Mei, R., “Cryogenic Two-Phase Flow During 
Chilldown,” Proceedings of the ASME HT-FED Conference, HT-FED2004-56555 Charlotte, NC, 
2004. 

Jackson, J., Liao, J., Klausner, J.F., Mei, R., “Transient Heat Transfer During Cryogenic 
Chilldown,” Proceedings of the ASME HT2005 Conference, HT2005-72145, 2005. 

Jackson, J., Klausner, J.F., and Mei, R., “Performance of Heat Transfer Correlations in 
Predicting the Heat Transfer During Cryogenic Chilldown,” Proceedings of the Sixth 
International Conference on Boiling Heat Transfer,” 2006. 


Students from Research 

1 . Chris Velat graduated with an M.S. degree. He is currently an Air Force Pilot 

2. Jelliffe Jackson graduated with a Ph.D. degree. He is currently a Professor at the University 
of Trinidad 

Nomeclature 

a cross-sectional area, m 2 

Bo Boiling number 

c, - c, constants in Kandlikar correlation 
Co Convection number 

e enhancement factor 

f„ fluid-dependent parameter in Kandlikar correlation 

f p pressure function in Mueller-Steinhagen correlation 

Fr Froude number 

G mass flux, kg/m 2 s 

m molecular weight, a. m.u. 

Pr Prandtl number, = c -^~ 

k 

Re Reynolds number 

r p surface roughness, m 
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5 suppression factor 

t temperature, °C 

x„ Martenilli parameter 

c p specific heat capacity, J/kgK 

d diameter, m 

/ friction factor 

g acceleration due to gravity, m/s 2 
h heat transfer coefficient, W/m 2 K 

h fg latent heat of vaporization, J/kg 

k thermal conductivity, W/mK 

n exponent in nucleate boiling correlation of Gorenflo 
p reduced pressure 

q heat flux, W/m 2 

r radial coordinate, m 

t time, s 

x quality 

z axial coordinate, m 

Greek Symbols 

<d annular flow modification factor in Muller-Steinhagen correlation 

8 liquid film thickness, m 

e scaling factor 

<!> azimuthal coordinate 

cp non-dimensional temperature 

p viscosity, Pa*s 

o dry dry angle of pipe perimeter, rad. 

p density, kg/m 3 

Subscripts 

o reference parameter 

00 far field 

dry dry region of pipe perimeter 

cb convective boiling 

fl fluid 

1 inner pipe surface 

known known 

/ liquid phase 

nb nucleate boiling 

out outer pipe surface 

pres present time instant 

prev previous time instant 

sat saturation condition 

tp two-phase 

v vapor phase 

w pipe wall 

wet wetted region of pipe perimeter 
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Superscripts 

‘ denotes non-dimensional variable 

~ denotes average quantity 
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APPENDIX A: Flow Boiling Heat Transfer Correlations 

Heat transfer coefficients from four well known nucleate flow boiling heat transfer correlations 
are used to compare against the experimentally determined average two-phase heat transfer 
coefficients defined in Eq. (8). Each correlation combines the contributions of the nucleate 
boiling and the convective mechanisms in different ways. 

The correlation of Gungor and Winterton [7] employs a superposition technique and accounts 
for the interaction of the two mechanisms through the use of an enhancement factor and a 
suppression factor. The enhancement factor accounts for the increased convective heat transfer 
that results from increased vapor fraction in the flow as it evaporates. The suppression factor 
reduces the nucleate boiling influence as more nucleation sites are suppressed with increasing 
flow velocity. The Gungor and Winterton correlation is as follows, 

h tp = Eh ,+Sh nb , (9) 

here \ is the two-phase heat transfer coefficient, h, is the convective heat transfer coefficient 
for the liquid flow, h nb is the nucleate boiling heat transfer coefficient, e is the enhancement 
factor and 5 is the suppression factor. The liquid convective heat transfer coefficient is 
calculated from the Dittus-Boelter correlation [11], 


h, =0.023 Re° 8 Pr, 0 ' 4 


V 

\ d iJ 


Re, is the liquid Reynolds number based on the liquid fraction flowing, 


A/ 


( 10 ) 


( 11 ) 


Here g is the total mass flux, * is the vapor quality, p is the viscosity, c p is the specific heat, 
k is the thermal conductivity, d t is the inner diameter of the pipe and the subscript / denotes 
the liquid phase. The nucleate boiling heat transfer coefficient is obtained from the Copper 
nucleate pool boiling correlation [12]: 

K =55^ 12 (-0.43431og^)- a55 M- 0 ^- 7 , (12) 

where p r is the reduce pressure and is computed as the ratio of the saturation pressure to the 
critical pressure of the fluid, m is the molecular weight of the fluid and q is the heat flux 
entering the system. The enhancement factor is a function of the Martinelli parameter, x tt and 
the Boiling number, Bo , 
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E = \ + 240005o +1.37 


where the Martinelli parameter and Boiling number are given as, 


(\~ 

9 f y 

Pv 

f X 

Pi 

V x J 

k F s j 

U j 


Here h fg is the latent heat of vaporization, p is density and the subscript v denotes the vapor 
phase. The proposed boiling suppression factor is, 


S' = (l + 1.15 xl0" 6 is 2 ReJ 17 ) _1 . 


In order to account for horizontal flows the enhancement factor is multiplied by a correction 
factor e 2 and the suppression factor is multiplied by a factor s 2 when the liquid Froude number, 
Fr, is such that Fr, < 0.05 , 


_ c, ( 0 . 1 — 2 Fr,) 


E 2 =Fr! 


S 2 ={Fr,)2, (18) 

where 

Fr i = —r~r- (19) 

Pi gd t 

The correlation of Kandlikar [8] accounts for the convective and nucleate boiling contributions to 
the heat transfer by introducing the Convection number, Co and the Boiling number, Bo . These 
quantities are used to enhance the single phase convective heat transfer coefficient; thus the 
correlation is often referred to as an enhancement model. The correlation takes the following 
form, 


K =(c i Co c - (25 Fr,f + C 3 Bo c < F fl )h„ 


where Co is, 
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The constants c,-c 5 are shown in Table 1 and are chosen based on the heat transfer regime 
for the given flow, and constant F fl depends on the fluid that is being utilized; ^ = 4.70 for 
nitrogen, other values may be found in [7]. 

Table 1. Constants used for the Kandlikar correlation 


Constant 

Convective 

region 

Nucleate 

boiling 

region 

c, 

1.1360 

0.6683 

c 2 

-0.9 

-0.2 

c 3 

667.2 

1058.0 

C 4 

0.7 

0.7 

c, 

0.3 

0.3 


C 5 = 0 for vertical tubes, and for horizontal tubes with Fr , > 0.04 . 


The correlation of Miiller-Steinhagen and Jamialahmadi [9] uses a superposition technique that 
accounts for the interaction of the convective and nucleate boiling mechanisms by employing an 
enhancement factor and a suppression factor. This formulation absorbs the enhancement factor 
into the computation of the convective heat transfer contribution, and the correlation takes the 
form, 

h tp =h l+ Sh nb . (22) 

This correlation employs the correlation of Petukov and Popov [13], for the liquid phase 
convective heat transfer coefficient, 


k i J A Rs„Pr, 


rf { 1+12 - 7 ^K J — ')) 


where the friction factor is given by Filonenko [14], 


/ = (1 .82 log Re, - 1 ,64)“ 2 . 


The two-phase Reynolds number is defined as 


Re, = G ( ! fH ffUs 


and the enhancement factor, e is given by, 
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f 


V 


E = 2.35 


— + 0.213 


v^« 

E = 1 otherwise 


J 


for > 0.1 

X„ 


For flow through annular tubes Eq. (23) is multiplied by a factor 


0 = 0.86 


f cl ^'°' 16 


\ d o J 


(26) 


(27) 


The suppression factor is computed as 


l + 2.53xlO“ 6 Re) n 17 " 

ip 


(28) 


In order to compute the nucleate boiling contribution to heat transfer, the pool boiling correlation 
developed by Gorenflo [15] is used, 



f \ 

q_ 

n 

\ Rp ) 

K*oJ 


^P° y 


(29) 


where r p is the surface roughness of the pipe. The pressure function, f p and the exponent n 
are calculated using the reduced pressure p , 


*0.27 

= 1 .73/? + 


6.1 + 


0.68 

1 * 2 

1 ~P 


(30) 


Equation (30) pertains to water and other low boiling point liquids and is used here; for the 
computation of f p for organic liquids reference [8] should be consulted. The exponent is 

calculated from 


n = 0.9 -0.3/“, (31) 

with a = 0.15 for water and other low boiling point liquids including nitrogen. The reference heat 
transfer coefficient h o , reference heat flux q 0 , and reference surface roughness R po for nitrogen 

W W 

are found in [16], and are q o =10000— , h a = 4380— — , and R =1x10 m . 

m m K 

The correlation that accounts for flow structure for horizontal flow is that of Wojtan et al. [6,10], It 
does so by first predicting the region of the pipe wall that is wetted and then computing the heat 
transfer for the wet and dry areas respectively. Only the heat transfer correlation is presented 
here. Determination of the flow structure is outlined in (10). The average two-phase heat 
transfer coefficient is predicted from 
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0 CyK + ( 2 X-O d , y )Ke, 


(32) 


where Q dry is the portion of the pipe cross section that is in contact with only the vapor; the 
determination of 9 dry is described in [6]. The convective heat transfer contribution by the vapor 
phase, 1 % , is computed using the Dittus-Boelter correlation, while the contribution of the liquid 
phase, h wet , is computed as 


(33) 

Here h cb is the convective boiling heat transfer coefficient and is based on a liquid film thickness, 
8 , length scale 


h = 

L wet 


[hj+{hj 


With 


h cb =0.0133Re®' 69 Pr / °' 4 ^- 



I ( d, V 2 A, 

V v 2 J ( 2 x-e dry ) 

otherwise. 


when 8 < — L 
2 


(34) 


(35) 


The nucleate boiling heat transfer coefficient is determined using the Copper nucleate pool 
boiling correlation [12], as is used in the work of Gungor and Winterton [7], and is modified by 
multiplying by a suppression factor of 0.8. 
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4. Numerical Investigation of Cryogenic Fluid Transport in Pipelines 
During Chilldown Process 
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Sub-Task 1: Numerical Investigation of Cryogenic Fluid Transport in Horizontal Pipelines 
for Low Mass-flux Regime during Chilldown Process 

Abstract 

A pseudo-steady model has been developed to predict the chilldown history of pipe wall 
temperature in the horizontal transport pipeline for cryogenic fluids. A new film boiling heat 
transfer model is developed by incorporating the stratified flow structure for cryogenic chilldown. 
A modified nucleate boiling heat transfer correlation for cryogenic chilldown process inside a 
horizontal pipe is proposed. The efficacy of the correlations is assessed by comparing the 
model predictions with measured values of wall temperature in several azimuthal positions in a 
well controlled experiment by Chung et al. (2004). The computed pipe wall temperature histories 
match well with the measured results. The present model captures important features of thermal 
interaction between the pipe wall and the cryogenic fluid, provides a simple and robust platform 
for predicting pipe wall chilldown history in long horizontal pipe at relatively low computational 
cost, and builds a foundation to incorporate the two-phase hydrodynamic interaction in the 
chilldown process. 

Nomenclature 

Bo Boiling number 

c solid heat capacity 

c p heat capacity 

D pipe diameter 

d thickness of pipe wall 

g gravity 

h heat transfer coefficient 

h poo i pool boiling heat transfer coefficient 

h conv convection boiling heat transfer coefficient 

h fg latent heat 

Ja Jacob number 

k thermal conductivity 

k eff effective thermal conductivity 

Nu Nusselt number 

p pressure 

R radius of pipe 

Rj and R 2 inner and outer radius of pipe 

Ra Rayleigh number 

Re Reynolds number 

Pc Peclet number 

Pr Prandtl number 
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5 

T 

T, 

T 2 

To 

t 

U 

u and v 
x andy 
Z 

z, r, and cp 
a 

Xtt 

S 

£ 

<P 

fl 

e 

a 


suppression factor 
temperature 
Leidenfrost temperature 

transition temperature between nucleate boiling to convection heat transfer 

room temperature 

time 

velocity 

vapor film velocity 
vapor film coordinates 
transformed coordinate 
cylindrical coordinates 
liquid volume fraction 
Martinelli number 
vapor film thickness 
emissivity 

azimuthal coordinate 
viscosity 

dimensionless temperature; azimuthal coordinate 
liquid surface tension; Stefan Boltzmann constant 


Subscripts 


0 characteristic value 

1 and o inner and outer pipe 

/ liquid 

v vapor 

w wall 

sat saturated 


Superscripts 

dimensionless variable 


I. Introduction 

The cryogenic chilldown is encountered in many applications but is of particular importance in 
cryogenic transportation pipelines. For example, in rockets or space shuttle launch facility, 
cryogenic liquids as fuel are filled from the storage tank to the internal fuel tanks of a space 
vehicle through a complex pipeline system. To avoid evaporated fuel entering the space vehicle, 
a cryogenic chilldown prior to the filling is required to reduce the pipe wall temperature to the 
saturation temperature of the cryogenic liquid. 

Cryogenic chilldown involves complicated hydrodynamic and thermal interactions among liquid, 
vapor, and solid pipe wall. There exist few basic experimental studies and modeling efforts for 
chilldown of cryogenic fluids. Studies on cryogenic chilldown started in 1960’s accompanying 
the development of rocket launching system. Early experimental studies were conducted by 
Burke et al. (ref. 1), Graham (ref. 2), Bronson et al. (ref. 3), Chi and Vetere (ref. 4), Steward 
(ref. 5) among other researches. Bronson et al. (ref. 3) studied the flow regime in a horizontal 
pipe during the chilldown by liquid hydrogen. The results revealed that the stratified flow is 
prevalent in the cryogenic chilldown. 

Flow regimes and heat transfer regimes in the horizontal pipe chilldown were also studied by 
Chi and Vetere (ref. 4). Information of flow regimes was deduced by studying the fluid 
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temperature and volume fraction during the chilldown. Several flow regimes were identified: 
single phase vapor, mist flow, slug flow, annular flow, bubbly flow, and single-phase liquid flow. 
Heat transfer regimes were identified as: single-phase vapor convection, film boiling, nucleate 
boiling, and single-phase liquid convection. 

Recently, Velat et al. (ref. 6) systematically studied cryogenic chilldown in a horizontal pipe. 
Their study included: visually recording the chilldown in a transparent Pyrex pipe, which is used 
to identify the flow regime and heat transfer regime; collecting temperature histories at different 
positions of wall in the chilldown; recording pressure drop along the pipe. Chung et al. (ref. 7) 
conducted a similar study on the nitrogen chilldown at relatively low mass flux and provided the 
data needed to assess various heat transfer coefficients in the present study. 

Burke et al. (ref. 1) developed a crude chilldown model based on the 1-D heat transfer through 
the pipe wall and the assumption of infinite heat transfer rate from the cryogenic fluid to the pipe 
wall. The effects of flow regimes on the heat transfer rate were neglected. Graham et al. (ref. 2) 
correlated heat transfer coefficient and pressure drop with the Martinelli number (ref. 8) based 
on their experimental data. Chi (ref. 9) developed a 1-D model for energy equations of the liquid 
and the wall, based on film boiling heat transfer between the wall and fluid. An empirical 
equation for predicting chilldown time and temperature was proposed. 

Steward (ref. 5) developed a homogeneous flow model for the chilldown. The model treated the 
cryogenic fluid as a homogeneous mixture. The continuity, momentum and energy equations of 
mixture were solved to obtain density, pressure and temperature of mixture. Various heat 
transfer regimes were considered: film boiling, nucleate boiling and single-phase convection 
heat transfer. Separate treatment of different heat transfer regimes resulted in a significant 
improvement in the prediction of the chilldown process. Homogeneous mixture model was also 
employed by Cross et al. (ref. 10) who obtained a correlation for the wall temperature in the 
chilldown with an oversimplified heat transfer model of the heat transfer between the wall and 
the fluid. 

Stratified flow regime, which is the prevalent flow regime in horizontal chilldown, was first 
studied by Chen and Banerjee (refs. 11 to 13). They developed a separated flow model for the 
simulating cooldown by a stratified flow in a hot horizontal pipe. Both phases were modeled 
using 1-D mass and momentum conservation equations. The wall temperature was computed 
using a 2-D transient heat conduction equation. Their prediction for wall temperature agreed 
well with their experimental results. Although a significant progress was made on handling the 
momentum equations, the heat transfer correlations employed were not as advanced. 
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Pipe wall 



Figure 1 . Schematic of chilldown heat transfer regimes in a horizontal pipe. 

Typical chilldown involves several heat transfer regimes as shown in figure 1. Near the liquid 
front is the film boiling. The knowledge of the heat transfer in the film boiling regime is relatively 
limited, because i) film boiling has not been the central interest in industrial applications; and ii) 
high temperature difference causes the difficulties in experimental investigations. For the film 
boiling on vertical surfaces, early work was reported by Bromley (ref. 14), Dougall and 
Rohsenow (ref. 15) and Laverty and Rohsenow (ref. 16). Film boiling in a horizontal cylinder 
was first studied by Bromley (ref. 17); and the Bromley correlation was widely used. Breen and 
Westwater (ref. 18) modified Bromley’s equation to account for very small tubes and large tubes. 
If the tube is larger than wavelength associated with the Taylor instability, the heat transfer 
correlation is reduced to Berenson’s correlation (ref. 19) for a horizontal surface. 

Empirical correlations for film boiling were proposed by Hendrick et al. (refs. 20 and 21), 
Ellerbrock et al. (ref. 22), von Glahn (ref. 23), Giarratano and Smith (ref. 24). These correlations 
relate a simple or modified Nusselt number ratio to the Martinelli parameter. Giarratano and 
Smith (ref. 24) gave detailed assessment on these correlations. All these correlations are for 
steady state cryogenic film boiling. They may not be suited for transient chilldown application. 

When the pipe chills down further, film boiling ceases and transient boiling occurs, followed by 
nucleate boiling. The heat transfer in the transient boiling is even more complicated and it is 
usually assumed that the boiling switches from film boiling to nucleate boiling right away. The 
position of film boiling transitioning to the nucleate boiling is often called rewetting front, 
because from that position the cold liquid starts touching the pipe wall. Usually the Leidenfrost 
temperature indicates the transition from the film boiling to the nucleate boiling. 

The study on the convection nucleate boiling is extensive. A general correlation for saturated 
boiling was introduced by Chen (ref. 25). Gungor and Winterton (ref. 26) modified Chen’s 
correlation and extend it to subcooled boiling. Enhancement and suppression factors for macro- 
convective heat transfer were introduced. Gunger and Winterton’s correlation can fit 
experimental data better than the modified Chen’s correlation (ref. 27) and Stephan and 
Auracher correlation (ref. 28). Kutateladze (ref. 29) and Steiner (ref. 30) also provided 
correlations for cryogenic fluids in pool boiling and forced convection boiling. Although they are 
not widely used, they should be more applicable for cryogenic fluids as the correlation was 
directly based on cryogenic conditions. 
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As the wall temperature drops further, the nucleate boiling is replaced by pure convection. The 
convection heat transfer can be modeled using Dittus-Boelter equation (ref. 31) for the fully 
developed turbulent pipe flow or corresponding laminar heat transfer equation for laminar pipe 
flow (ref. 31). The condition in which nucleate boiling switches to single-phase convection heat 
transfer is that all nucleate sites are suppressed. 

Although two-fluid model can describe the fluid dynamics aspect of the chilldown process, it 
suffers from computational instability for moderately values of slip velocity between two phases 
which limits its application. To gain fundamental insight into the thermal interaction between the 
wall and the cryogenic fluid and to be able to rapidly predict chilldown in a long pipe, an 
alternative pseudo-steady model is developed. In this model, liquid wave front speed is 
assumed to be constant and is the same as the bulk liquid speed (ref. 32). It is also assumed 
that steady state thermal fields for both the liquid and the solid exist in a reference frame that is 
moving with the wave front. The heat transfer between the fluid and the wall is modeled using 
different heat transfer correlations depending on the operating heat transfer regime at a given 
location. Various improvements on the correlations are introduced, including the development of 
a new film boiling heat transfer coefficient. The governing equation for the solid thermal field 
becomes a parabolic equation that can be efficiently solved. It must be emphasized that a great 
advantage of the pseudo-steady model is that one can assess the efficacy of the film boiling 
model independently from that of the nucleate boiling model since the down stream information 
in the nucleate boiling regime cannot affect the temperature in the film boiling regime. In another 
word, even if the nucleate boiling heat transfer coefficient is inadequate, the film boiling heat 
transfer coefficient can still be assessed in the film boiling regime by comparing with the 
measure temperature for the right period of time. After the satisfactory performance is achieved 
for the film boiling regime, the nucleate boiling heat transfer model can be subsequently 
assessed. In the Results section, those detailed assessments of the heat transfer coefficients 
are provided by comparing the computed temperature variations with the experimental 
measurements of Chung et al. (ref. 7). Satisfactory results are obtained. 

II. Formulation 

In pseudo-steady chilldown model, it is assumed that both the liquid and its wave front moves at 
a constant speed U which is taken from estimated experimental condition in the present model. 
Thus the main emphasis of the present study is on the modeling of the heat transfer coefficients 
in different heat transfer regimes and the computation of the thermal field within the solid pipe. 

11.1. Solid Heat Transfer 


The thermal field inside the solid wall is governed by the 3-D unsteady energy equation: 


8T_ _d_, 
dt dz 


k— 

dz 


1 d 

H 

r dr 


' , dT\ 1 d (kdT) 
V dr ) r d(p\r d(p J 


( 1 ) 


Since the wave front speed U is assumed to be a constant, it can be expected that when the 
front is reasonably far from the entrance region of the pipe, the thermal field in the solid is in a 
steady state when it is viewed in the reference frame that moves with the wave front. Thus, the 
following coordinate transformation is introduced, 

Z = z + Ut (2) 
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Equation (1) is then transformed to: 


tt 8T d ( 7 dT 
pcU — = — k — 
dZ dZ{ dZ 


1 d 
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r dr 


, dT} 1 d (kdT\ 
dr ) r d(p \ r d(p , 


(3) 


For further simplification, the following dimensionless parameters are introduced, 
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T-T 
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T -T 
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r , c , k 

— c = — and k = — 

d Cq kq 


(4) 


where k 0 is the characteristic thermal conductivity, and c 0 is the characteristic heat capacity. 


Equation (3) can be normalized as 


Pc * 



d 

( k , S 0) 

1 d 

\ 

r k — 

dZ' 

{ SZ'J 

r' dr' 

l dr'J 


1 d 

r' dcp 


f iP_de\ 

r' (p 


(5) 


pcJJd 

where Pc = — - — is the Peclet number. 

K 

Under typical operating condition for the cryogenic chill-down process, Pc ~ 0(1 0 2 ) - 0(1 0 3 ). 
The first term on the RHS of Eq. (5) is quite small compared with the rest of terms and can be 
neglected. Eq. (5) becomes 


Pc*c' 


dQ_ 

dZ’ 


j__a_ 

r' dr' 


r , h ,de' 
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dr' 


1 _AS—— 

r' dtp yr' q> J’ 


( 6 ) 


which is parabolic. Hence in the Z'-direction, only one condition is needed. In the(p-direction, 
periodic boundary conditions are used. On inner and outer surfaces of the pipe wall, where r = 
R 1 and R 2 , proper boundary conditions for the temperature are required. 



Figure 2. Coordinate systems: laboratory frame is denoted by z; moving frame is denoted by Z. 
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For convenience, Z' = 0 at the liquid wave front. In the region of Z' <0 , the inner wall is 
exposed to the pure vapor. Even though there may be some liquid droplets in the vapor that can 
cause evaporative cooling when the droplets deposit on the wall and the cold flowing vapor can 
absorb some heat from the wall, the heat flux due to these two mechanisms is quit small 
compared with the heat transfer between liquid and solid wall in the region of Z' > 0 . Hence, 
heat transfer for Z' > 0 is neglected and it is assumed that 9 = 1 at Z' = 0. The computation 
start from the Z' = 0 to Z' — > oo , until a steady state solution in the Z' -direction is reached. An 
implicit scheme in the Z'-direction is employed to solve eq. (6). 

11.2. Liquid and Vapor Flow 

The two-phase flow is assumed to be stratified as was observed in (ref. 7). Both liquid and 
vapor phases are assumed to be at the saturated state. Liquid volume fraction determining 
which part of wall is in contact with the liquid or the vapor, is specified at every cross-section 
along the Z'-direction based on experimental information. The visual studies (refs. 6 and 7) 
show that the liquid volume fraction increases gradually, rather than abruptly, near the liquid 
wave front and becomes almost constant during most of the chilldown. Hence, the following 
liquid volume fraction a as a function of time is used in the computation of the solid-fluid heat 
transfer coefficient, 


a = a n sin 

U 2j 

t<t 0 


(7) 

a - a 0 

IV 

O 



where t 0 is characteristic chilldown time, andcr 0 is characteristic liquid volume fraction. Here 

the time when the nucleate boiling is almost suppressed and slope of the wall temperature 
profile becomes flat is set as characteristic chilldown time. 

The vapor phase velocity is assumed to be a constant. However, it was not directly measured. It 
is computationally determined by trial-and-error by fitting the computed and measured wall 
temperature variations for numerous positions. 

11.2. 1. Heat Transfer Between Cryogenic Fluid and Solid Pipe Wall 

During the chilldown, the fluid in contact with the pipe wall is either liquid or vapor. The 
mechanisms of heat transfer between the liquid and the wall and between the vapor and the 
wall are different. Based on experimental measurements and theoretical analysis, liquid-solid 
heat transfer accounts for a majority of the total heat transfer. However, the prediction for heat 
transfer is much more complicated than the heat transfer between the vapor and the wall. The 
heat transfer between the liquid and the wall is discussed first. 

11.2.1.1. Heat Transfer Between Liquid and Solid wall 

Heat transfer between liquid and solid wall includes film boiling, nucleate boiling, and single- 
phase convection heat transfer. The transition from one type of boiling to another depends on 
many parameters, such as wall temperature, wall heat flux, and various properties of fluid. For 
simplicity, a fixed temperature approach is adopted to determine the transition points. That is, if 
the wall temperature is higher than the Leidenfrost temperature, film boiling dominates. If the 
wall temperature is between Leidenfrost temperature and a transition temperature, T 2 , nucleate 
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boiling takes place. Here, the transient boiling regime between film boiling and nucleate boiling 
is neglected. The reason is mainly the difficulty associated with the determination of the 
conditions for the transient boiling to occur and the lack of a reliable transient heat transfer 
correlation. If the wall temperature is below the transition temperature T 2 , pure convection heat 
transfer dominates. The values of Leidenfrost temperature and transition temperature are 
determined by matching model prediction with the experimental results. 

11.2.1.1.1. Film boiling heat transfer. 

Because of the nature of high wall superheat in chilldown, the film boiling plays a major role in 
chilldown in terms of the time span and in terms of the total amount of heat removed from the 
wall. Currently there exists no specific film boiling correlation for chilldown applications with such 
high superheat. If one uses conventional film boiling correlations, necessary modifications for 
cryogenic application must be made for chilldown. 

One of cryogenic film boiling heat transfer correlations was provided by Giarratano and Smith 
(ref. 24) 


f Nu ' 

\ NU calcJ 


Bo~ 0A = f(Ztt) 


( 8 ) 


where Nu calc is the Nusselt number for the forced convection heat transfer. In this correlation, 

the heat transfer coefficient is the averaged value for the whole cross section. Similar 
correlations for cryogenic film boiling also exist in literature. The correlations were obtained from 
measurements conducted under steady states. The problem with the use of this steady state 
film boiling correlation is that it does not take into account the change of flow regime as 
encountered in the chilldown. For example, for the same quality, the heat transfer rate in 
annular flow is much different from that in stratified flow, while those empirical correlations 
cannot take such difference into account. 

Furthermore, in this study, local heat transfer coefficient is needed in order to incorporate the 
thermal interaction with the pipe wall. Since the two-phase flow regime information is available 
through visualization in the present study, the modeling effort needs to take into account the 
knowledge of the flow regimes. 

There are several correlations for the film boiling based on the analysis of vapor film boundary 
layer and stability of the thin vapor film, such as Bromley’s correlation (ref. 17) and Breen and 
Westerwater’s correlation (ref. 18) for film boiling on the outer surface of a hot tube, Frederking 
and Clark’s (ref. 33) and Carey’s (ref. 34) correlations for film boiling on the surface of a sphere. 
However, none of these was obtained for cryogenic fluids or for film boiling on the inner surface 
of a pipe. 

A new correlation for the film boiling in cryogenic chilldown inside a tube is presented here. The 
schematic of film boiling inside a pipe is shown in figure 3 with a cross-sectional view. Bulk 
liquid is near the bottom of the pipe. Beneath the liquid is a thin vapor film. Due to buoyancy 
force, the vapor in the film flows upward along the azimuthal direction. Heat is transferred 
through the thin vapor film from the solid to the liquid. Reliable heat transfer correlation for film 
boiling in pipes or tubes requires the knowledge of the thin vapor film thickness which can be 
obtained by solving the film layer continuity, momentum, and energy equations. 
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To simplify the analysis for the vapor film heat transfer, it is assumed the liquid velocity in the 
azimuthal direction is zero and the vapor flow in the direction perpendicular to the cross-section 
is negligible. It is further assumed that the vapor film thickness is small compared with the pipe 
radius and vapor flow is in steady state, incompressible and laminar. The laminar flow 
assumption can be confirmed post priori as the Reynolds number, Re, based on the film velocity 

and film thickness is typically of o(l0° ~ 10 2 ). In terms of the x- and y-coordinates and (u, v) 
velocity components shown in figure 3, the governing equations for the vapor flow are similar to 
boundary-layer equations: 


du dv 
dx dy 


( 9 ) 


du du 1 dp d 2 u . . 

u — + v — = + v v — -- esin#, 

dx dy p v dx dy 2 


( 10 ) 


dr dr dr 

u h v — = a v — - 

dx dy dy 2 


( 11 ) 


where subscribe v represents properties of vapor. 

Because the length scale at x direction is much larger than the length scale at y direction, the v 
velocity can be neglected. Furthermore, convection term is assumed small and is neglected, 
and the momentum equation is simplified to 


1 dp 
P v dx 



- g sin <9 . 


( 12 ) 


The vapor pressure can be evaluated by considering the hydrostatic pressure from liquid core 
as: 
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p = p 0 + p,gR cos — -cos^ 
1 R / 


where $ is angular position where the film merges with the vapor core. The momentum equation 
becomes 


(a-a)_ . , .. 3 


■ gsin UJ +,/ -^ =0 


The vapor velocity boundary condition is u = 0 at y = 0 and u = u, = 0 at y = 8 . The vapor 
velocity profile is: 

< 15 > 

2 v v p v \R) 


The mean u velocity is defined as 


- i r* j ( Pi ~ Pv)<5 2 g ■ f x 
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The energy and mass balance on the vapor film requires that 


l By 


= dm = p v d(uS) 


Neglecting the convection, the vapor energy equation is: 


The following linear temperature profile is thus obtained, 


T-T„ t y 
T,~T m S 


Substituting the temperature and velocity profiles into eq. (17) yields 
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Equation (20) has analytical solution: 
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where Ja is Jacob number and Ra is Raleigh number: 
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describes the geometric dependence of vapor film thickness. 
The mean velocity u as a function of 6 is thus 


u = 


C T w - T sat\pl ~Pv)g R 
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( 25 ) 
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Figure 4. Numerical solution of the vapor thickness and velocity influence functions. 


Curves for F(0 ) andF(#) 2 sin# based on numerical integration are shown in figure 4. The 


n 


vapor film thickness has a minimum at # = 0 and is nearly constant for 6 < — . It rapidly grows 


after # > — . The singularity at the top of tube when # — » n is of no practical significance since 

the film will merge with the vapor core at the vapor-liquid interface. The vapor velocity is 
controlled by F(#) 2 sin# which is zero at the bottom of the pipe and increases almost linearly in 
the lower part of the tube where the vapor film thickness does not change substantially. In the 
upper part of the tube, due to the increase in the vapor film thickness, the vapor velocity 
gradually drops back to zero at the top of the tube. Thus a maximum velocity may exist in the 
upper part of the tube. 


The local film boiling heat transfer coefficient is easily obtained from the linear temperature 
profile. It is 


h = 


K 
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0.6389 
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11.2.1.1.2. Nucleate Boiling and Convection Heat Transfer 


(26) 
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For the nucleate flow boiling heat transfer, Gungor and Winterton’s correlation (ref. 26) is widely 
used due to its relatively better accuracy in predicting heat transfer coefficients. However, a 
closer examination on this correlation shows that it is based mainly on the following parameters: 
Pr, Re, and quality x. Similar to the development of conventional film boiling correlations, these 
parameters all reflect overall properties of the flow in the pipe and are not directly related to flow 
regimes. Thus it cannot be used to predict the local heat transfer coefficients in chilldown. 
Chen’s (ref. 25) flow boiling correlation is based on separating the heat transfer to micro- and 
macroconvection heat transfer. Micro-convection heat transfer represents the contribution from 
boiling heat transfer, and macro-convection represents the contribution from the forced 
convection heat transfer. However his correlation fits best for annular flow because of the 
derivation of suppression factor. In stratified flow regime, which is common in cryogenic 
chilldown, Chen’s correlation may not be applicable. 

Several correlations have been tried in this study, including Gungor and Winterton’s correlations 
(ref. 26), Chen’s correlation (ref. 25), and Kutateladze’s correlations (ref. 29). None gives a 
satisfactory heat transfer rate that is need to match the experimentally measured temperature 
histories in (ref. 7) in the nucleate boiling regime. Among them, Kutateladze’s correlation gives 
more reasonable result. In this correlation, the total heat transfer coefficient h is 


h = h + h 


pool 


(27) 


where h conv is given by Dittus-Boelter equation for fully developed pipe flow, 


A„„,=0.023*Re?- 8 Pr,°- 4 *t,/Z) 


(28) 


And h pool is pool nucleate boiling heat transfer coefficient, 


h 


pool 


0.487 *1CT 10 * 


/ 1 .282 1 .750 / V - 5 

k iPi P { Cp), 

77 „ \l -5 „ 0.906 ,, 0.626 

\h fg p v ) o p, 


AT 15 


(29) 


In which AT is wall superheat. 

Kutateladze’s correlation (ref. 29) was proposed without considering the effect of nucleate site 
suppression. This obviously leads to an overestimation of the nucleate boiling heat transfer rate. 
Hence a modified version of Kutateladze’s correlations (28) is used, 

h =Konv+ S * h pool ( 30 ) 


where S is suppression factor. 

When AT drops to a certain range all the nucleate sites are suppressed. The heat transfer is 
dominated by single phase forced convection. The heat transfer coefficient can then be 
predicated using Dittus-Boelter equation, eq. (27), if flow is turbulent, or eq. (30), if flow is 
laminar. 


K onv = 4.36* k,/D 


(31) 
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11.2.1.2. Heat Transfer between Vapor and Solid Wall 

The heat transfer between the vapor and the wall can be estimated by treating the flow as a fully 
developed convection flow, neglecting the liquid droplets that are entrapped in the vapor. The 
heat transfer coefficient of vapor forced convection flow is 

h v = 0.023 * Re“ 8 Pr v ° 4 *k v /D, (turbulent flow) (32) 


or 


h v =4.36 *k v /D, (laminar flow) (33) 

11.2.2. Heat Transfer between Solid Wall and Environment 

For a cryogenic flow facility, although serious insulation is applied, the heat leakage to 
environment is still considerable due to the large temperature difference between the cryogenic 
fluid and the environment. It is necessary to evaluate the heat leakage from the inner pipe to 
environment in order to make realistic assessment of the model prediction with the experimental 
results (ref. 7). 

Vacuum insulation chamber between the inner and outer pipes is used in cryogenic transport 
pipe (ref. 7), as shown in figure 5. Radiation heat transfer exists between the inner and outer 
pipe. Furthermore, the space between the inner and outer pipe is not an absolute vacuum. 
There is residual air that causes free convection between the inner and outer pipes driven by 
the temperature difference of the inner and outer pipe. 

The radiation between the inner pipe and outer pipe becomes significant when the inner pipe is 
chilled down. The heat transfer coefficient is proportional to the difference of the fourth power of 
wall temperatures. Exact evaluation of the radiation heat transfer between the inner and outer 
pipe is a difficult task. Hence a simplified model based on the overall radiation heat transfer 
between long concentric cylinders with constant temperature at inner pipe and outer pipe 
(ref. 31) is used to evaluate the heat transfer rate at every axial location of the pipe. It is not 
quantitatively correct, but can provide reasonable estimate for the magnitude of the radiation 
heat transfer between pipes through the vacuum. The local radiation heat transfer rate per unit 
area on the surface of inner pipe rad q' radm is 


Qrad ~ 


~K) 


1 1 -i 

— + 

£; £„ 


r r ' 


\ r oJ 


(34) 


where the o is Stefan Boltzmann constant, (n, si) and (r 0 , s 0 ) are the radius and emissivity of 
inner pipe and outer pipe, respectively, and T wa ii is the local inner wall temperature, T 0 is the 
room temperature that is assumed constant in the entire outer pipe. Here the emissivity is also 
assumed to be constant during the entire chilldown process. 

For the free convection heat transfer due to the residual air in the vacuum chamber between the 
inner pipe and outer pipe, Raithby and Hollands’ correlation (ref. 35) is used for the heat transfer 
rate. The average heat transfer rate per unit length of the cylinder is 
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(T,-T 0 ) 


(35) 
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eff 
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where T is assumed constant on the inner and outer wall along the azimuthal directions, and eff 
k is the effective thermal conductivity. 



Figure 5. Schematic of vacuum insulation chamber. 

III. Results and Discussions 

In experiment by Chung et al. (ref. 7), liquid nitrogen was used. The flow regime is revealed to 
be stratified flow by visual observation, as shown in figure 7, and wall temperature history in 
several azimuthal positions is measured by thermal couples and recorded on a computer. 

111.1. Experiment of Chung et al. 

In experiment of Chung et al. (ref. 7), a concentric pipe test section (fig. 6) was used. The 
chamber between the inner and outer pipe is vacuumed, but about 20 percent air remained. The 
inner diameter (i.d.) and outer diameter (o.d.) of inner pipe are 11.1 and 15.9 mm, and i.d. and 
o.d. of outer pipe are 95.3 and 101.6 mm, respectively. Numerous thermal couples were placed 
at different locations of the inner pipe. Some were embedded very closely to the inner surface of 
the inner pipe while others measure the outside wall temperature of inner pipe. Experiments 
were carried out at room temperature and atmosphere pressure. Liquid nitrogen flows from a 
reservoir to the test section driven by gravity. As liquid nitrogen flows through the pipe, it 
evaporates and chills the pipe. Some of the typical visual results are shown in figure 7. The 
measured average liquid nitrogen velocity is U ~ 5 cm/s. Vapor velocity is not measured in this 
experiment. In this study, it is determined through trial-and-error by fitting the computed and 
measured temperature histories. The characteristic liquid volume fraction is 0.3 from the 
recorded video images. The characteristic time used in this computation is t 0 =100s. The 

Leidenfrost temperature for the nitrogen is around 180 K; hence the temperature when the film 
boiling ends and nucleates boiling starts is set at 180 K. The transition temperature at which the 
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nucleate flow boiling switches to purely convection heat transfer is 140 K based on experimental 
results. The material of the inner pipe and outer pipe used in the experiment of (ref. 7) are Pyrex 
glass with emissivity of 0.82 at room temperature. 


Slordord To 



Figure 6. Schematic of Chung et al.’s cryogenic twophase flow test apparatus (ref. 7). 



Figure 7. Experimental visual observation of Chung et al.’s cryogenic two-phase flow test (ref. 7). 
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III. 2. Comparison of pipe wall temperature 

In the computation, there are 40 grids along the radial direction and 40 grids along azimuthal 
direction for the inner pipe (fig. 8). The results of temperature profile at 40-40 grids and higher 
grids resolution show that 40-40 grids are sufficient. Figures 9 and 10 compare the measured 
and computed wall temperature as a function of time at positions 11, 12, 14, and 15 shown in 
figure 8. For the modified Kutateladze’s correlation a suppression factor 0.01 is used. Vapor 
velocity is 0.5 m/s based on the best fit. 

The overall temperature histories agree well in the film boiling stage. Thus, the film boiling heat 
transfer coefficient based on the first principle and incorporated the flow structure gives very 
reliable prediction. It must also be noted that the value of the Leidenfrost temperature does not 
affect the computed temperature history prior to that transition point since the governing 
equation is parabolic. The good agreement before the Leidenfrost temperature is reached is 
entirely due to the superior performance of the new film boiling heat transfer coefficient. 

During the stage of the rapidly decreasing wall temperature after the Leidenfrost temperature, 
the computed wall temperature drops slightly faster than the measured value. The rapid 
decrease in the wall temperature is due to initiation of nucleate boiling which is more efficient for 
heat transfer than the film boiling. Reasonable agreement between the computed and measured 
histories in this nucleate boiling regime is due to: i) the good agreement already achieved in the 
film boiling stage; ii) valid choice for the Leidenfrost temperature that switches the heat transfer 
regime correctly; and iii) appropriate modification of Kutateladze’s correlations. 

In the final stage of chilldown, the wall temperature decreases slowly, and the computed wall 
temperature shows the same trend as the measured one but tends to be a little lower. Figure 1 1 
shows the temperature distribution of a given cross-section at different times of chilldown. 
Because the upper part of pipe wall is exposed to nitrogen vapor, cooling effect is much 
reduced. 



Figure 8. The computational grids arrangement and positions of thermal couples. 
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Figure 9. Comparison of wall temperature of transducer 12 and 15, which is at the bottom of 
pipe, between the experiment and computation. For comparison purpose, the solution in the 

Z- direction is converted to time t. 



Figure 10. Comparison of wall temperature of transducer 1 1 and 14, which are at the upper part 
of the pipe, between the experiment and computation. For comparison purpose, the solution in 

the Z-direction is converted to time t. 
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t=Os 


t=100s 



t=50s 


I 


T 

262.027 
258.192 
254.356 
250.521 
246.685 
242.85 
239.014 
235.179 
231 .343 
227.507 
223.672 
219.836 
216.001 
212.165 
208.33 


1 229.257 
222.665 
216.073 
209.481 
202.89 
196.298 
189.706 
183.114 
176.522 
169.93 
163.338 
156.746 
150.155 
143.563 
136.971 


t=300s 


0.006 

0.004 

0.002 

>■ 0 


- 0.002 


-0.004 


-0.006 



1 149.383 
147.146 
144.91 
142.673 
140.437 
138.2 
135.964 
133.727 
131.491 
129.254 
127.018 
124.781 
122.545 
120.308 
118.072 
115.835 
113.599 

109.126 

106.889 




Figure 1 1 . Wall temperature distribution across the section at t = 0, 50, 100, and 300 sec. 

III. 3. Remarks 

In figure 9, the wall temperature based on the film boiling correlation of Giarratano and Smith 
(ref. 24) is also shown. Apparently, the correlation of Giarratano and Smith (ref. 24) gives a very 
low heat transfer rate so that the wall temperature remains high. This comparison confirms our 
earlier argument that correlations based on the overall flow parameter, such as quality and 
averaged Reynolds number, are not applicable for the simulation of unsteady chilldown. 

The nucleate flow boiling correlations of Gungor and Winterton (ref. 26), Chen (ref. 25), 
Kutateladze’s correlations (ref. 29) are also compared in this study. Gungor and Winterton’s 
correlation fails to give a converged heat transfer rate during the transition from the film boiling 
to the nucleate boiling. Chen’s correlation overestimates the heat transfer rate, and cause 
unrealistically large temperature drop on the wall, which results in the halt of the computation. 
Only Kutateladze’s correlation gives a reasonable heat transfer rate. However, the temperature 
drop near the bottom of the pipe is still faster than the measured one as shown in figure 9. This 
may be due to the fact that most of nucleate boiling correlations were obtained from 
experiments of low wall superheat. However, in this cryogenic chilldown, wall superheat is much 
higher than that in the normal nucleate boiling experiments. Another reason is that the original 
Kutateladze’s correlation does not include suppression factor. This leads to overestimating heat 
transfer coefficient. The modified Kutateladze’s correlation with suppression factor gives 
reasonable chilldown result in figure 9. 

Further examination of figures 9 and 10 indicates that although we have considered the heat 
leak from the outer wall to the inner wall through radiation and free convection, the computed 
temperature is still lower than the measured temperature during the final stage of chilldown, 
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where the heat transfer rate between fluid and wall is low due to lower wall superheat and 
possibly heat leak. The temperature difference between the computed and measured values at 
position 12 and position 15 suggests that there may be additional heat leak, which affects the 
measurements but is not taken in account in the present modeling. 

In this study, the pseudo-steady chilldown model is developed to predict the chilldown process 
in a horizontal pipe in stratified flow regime. This model can also be extended to describe 
annular flow chilldown in the horizontal or the vertical pipe with minor changes on the boundary 
condition for solid temperature. It can also be extended to study the chilldown in the slug flow as 
long as we specify the contact period between the solid and the liquid or the vapor. The 
disadvantage of the current pseudosteady chilldown model is that the fluid interaction inside the 
pipe is largely neglected and both vapor and liquid velocities are assumed to be constant. 
Compared with a more complete model that incorporates the two-fluid model, the present 
pseudo-steady chilldown model requires more experimental measurements as inputs. However, 
the pseudo-steady chilldown model is computationally more robust and efficient for predicting 
chilldown process. It provides overall reasonable results for the solid wall temperature. While a 
more complete model for chilldown process that incorporates the mass, momentum, and energy 
equations of vapor and liquid is being developed to reduce the dependence of the experimental 
inputs for the liquid velocity and trial-and-error for the vapor velocity, the present study has 
revealed useful insights into the key elements of two-phase heat transfer encountered in the 
chilldown process which have been largely ignored. It also laid the necessary modeling 
foundation for the incorporation of the two-fluid model. 

IV. Conclusions 

A pseudo-steady chilldown computational model has been developed to understand the heat 
transfer mechanisms of cryogenic chilldown and predict the chilldown wall temperature history 
in a horizontal pipeline. The model assumes the constant speed of the moving liquid wave front, 
and steady state of the thermal field in the solid in a reference frame that moves with the liquid 
wave front, and saturated temperature of both liquid and vapor so that the 3-D unsteady 
problem can be transformed to a 2-D, parabolic problem. A new film boiling heat transfer 
coefficient in the cryogenic chilldown condition is developed using the first principle and 
incorporating the stratified flow structure. The existing nucleate boiling heat transfer correlations 
may not work well under cryogenic condition. A modified Kutateladze’s correlations with 
suppression factor adequately describes heat transfer coefficient. With the new and modified 
heat transfer correlations, the pipe wall temperature history based on the pseudo-steady 
chilldown model matches well with the experimental results by Chung et al. (ref. 7) for almost 
the entire chilldown process. The pseudo-steady chilldown model has captured the important 
features of thermal interaction between the pipe wall and the cryogenic fluid. 
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Sub-Task 2: A Study on Numerical Instability of Inviscid Two-Fluid Model near 
lll-posedness Condition 


Abstract 

Two-fluid model is widely used in studying gas-liquid flow inside pipeline, because it can 
qualitatively predict the flow field at low computational cost. However, the two-fluid model is 
quite unstable when it reaches or is near ill-posedness. In this study a pressure correction 
algorithm for inviscid two-fluid model is carefully designed to improve stability near ill-posedness. 
Through von Neumann stability analysis, the wave growth rates by pressure correction scheme 
with first order upwind, second order upwind, and central difference dicretization scheme are 
studied. The von Neumann analysis shows that central difference scheme is more accurate and 
stable than other two schemes. Compared with first order upwind scheme, second order upwind 
scheme, which is supposed to be more accurate and widely used by CFD community, is much 
unstable for low frequent harmonics and inaccurate for high frequency harmonics in two-fluid 
model. Von Neumann stability analysis also shows the instability by ill-posedness of two-fluid 
model is significantly different from the instability of the discretized two-fluid model. Furthermore, 
comparisons between actual wave growth and predicted wave growth by von Neumann stability 
analysis are taken. The comparison shows the accuracy of the new pressure correction scheme 
and success of von Neumann stability analysis in two-fluid model. Accordingly, the relation 
between the ill-posedness of two-fluid model and numerical stability of algorithm for inviscid two- 
fluid model is elucidated. 


Nomenclature 

c wave speed 

E common amplitude factor 

G amplification factor 

g gravity 

k wavenumber 

H hydraulic depth 

p pressure 

N grids number 

t time 

u velocity 

x space coordinate 

a volume fraction 

p angle of inclination from the horizontal 

s amplitude factor 

cp phase angle 

X characteristic root 

p denisty 

Subscripts 

e east face of control volume 

g gas 

/ liquid 

i interface, grid index 

P center of main control volume 

w west face of control volume 
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I. Introduction 


Gas-liquid flow inside a pipeline is prevalent in the handling and transportation of fluid. A reliable 
flow model is essential to the prediction of flow field inside the pipeline. To fully simulate the 
system, Navier-Stokes equations in three-dimensions are required. However, it is not feasible 
for simulating flow in a long pipe by today’s computer. To reduce the computational cost in 
simulating the flow in a long pipe and obtain basic and essential flow properties, such as gas 
volume fraction, liquid and gas velocity, pressure, one-dimensional model is necessary. Two- 
fluid model as one of one-dimensional model is recognized as a realistic tool to simulate the 
gas-liquid flow inside the pipeline. Two-fluid model is also called separated flow model, which 
consists of two sets of conservation equation for mass, momentum and energy of gas phase 
and liquid phase. It was proposed by Willis (1969), and further refined in Ishii (1975). 

Although two-fluid model shows its success in simulating the two-phase flow in pipeline, the 
two-fluid model suffers ill-posedness problem that is, when relative velocity between liquid and 
gas exceeds a critical value, the governing equations do not possess real characteristic 
(Gidaspow, 1974; Jones and Prosperettii, 1985; Song and Ishii, 2000). This ill-posedness 
condition suggests that the results of two-fluid model at that condition do not reflect real flow 
situation inside the pipe. The two-fluid model only gives meaningful results when the relative 
velocity between gas and liquid phase is less than a criterion, which is determined by gravity, 
liquid level, surface tension, and other conditions. However, the criterion is coincident with the 
stability criterion of inviscid Kelvin-Helmholtz instability (IKH) analysis (Issa and Kempf, 2002). 
Because the instability of IKH analysis results in the flow regime transition from stratified flow to 
slug flow or annular flow (Barnea and Taitel, 1994a), ill-posedness of two-fluid model also 
implies the flow regime transition (Brauner and Maron, 1992). 

The computational methods to solve governing equations of two-fluid model were researched by 
many investigators. In this study, a further assumption on two-fluid model is both liquid and gas 
phases are incompressible, which is reasonable because most of stratified flows are at relative 
low speed compared with sound speed. However, if the fluids are assumed incompressible, the 
governing equations of two-fluid model are not pure hyperbolic. Therefore hyperbolic type solver 
cannot be used on incompressible two-fluid model. To solve the incompressible two-fluid model, 
one approach is to simplify the governing equations to two governing equations about liquid 
phase volume fraction and liquid velocity, but transient terms in gas mass and momentum 
equations are neglected (Chan and Banerjee, 1981; Barnea and Taitel, 1994b). The governing 
equations of simplified incompressible two-fluid model are hyperbolic, and can be solved by 
characteristic method or finite difference method. Another more effective method of solving 
incompressible two-fluid model is to use pressure correction scheme (Patanka 1980). Issa and 
Kempf (2000) successfully applied the pressure correction scheme in two-fluid model and 
successfully simulated stratified flow and slug flow inside a pipe. 

When two-fluid model becomes ill-posed, the solution of the model becomes unstable. The 
discretized model should demonstrate the instability of two-fluid model, but numerical instability 
may not be identical to the instability caused by the ill-posedness. Inappropriate numerical 
schemes may trigger the numerical instability earlier than the actual one. Lyczkowski et al. 
(1978) used von Neumann stability analysis to study a compressible two-fluid model with their 
numerical scheme. Their study shows numerical instability and ill-posedness may not be 
identical. However, their study is incomplete, lacking gravity term and analysis on different 
discretization schemes. Stewart (1979), Ohkawa and Tomiyama (1995) tried to illustrate the 
numerical stability of incompressible two-fluid model with a simplified model equation as a 
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replacement. Their study shows that higher order upwind schemes yield more unstable 
numerical solution than the first order upwind scheme. 

A complete study on the numerical stability of two-fluid model and the relation between 
numerical instability and ill-posedness are desired. 

In this study, a pressure correction scheme designed to increase the stability of numerical 
scheme when flow is near ill-posedness will be present. Von Neumann stability analysis is 
employed to study the stability of discretized two-fluid model with different interpolation schemes. 
The performances of different schemes at various flow configurations will be compared. The 
results of von Neumann stability analysis can be used to validate the computer code. 
Comparisons will be taken at various flow configurations. Furthermore, test of wave propagation 
inside pipe shows the scheme has the ability to handle the real flow in the pipe. 

Ill-posedness of two-fluid model is the major concern in this study. Since ill-posedness and IKH 
stability are not related to the viscosity of two-fluid model. To simplify the analysis on the two- 
fluid model and concentrate our major concern, only inviscid two-fluid model are considered. 
Furthermore all of results of study on the inviscid two-fluid model are also applicable to the 
viscous two-fluid model. 

II. Numerical Method 

11.1. Governing Equations 

The basis of the two-fluid model is a set of one-dimensional conservation equations for the 
balance of mass, momentum and energy for each phase. The one-dimensional conservation 
equations are obtained by integrating the flow properties over the cross-sectional area of the 
flow. 



Pipe cross section 



Fig.1 Schematic of two-fluid model for pipe flow. 


Because the ill-posedness originates from the hydrodynamic instability of the two-fluid model, 
only continuity and momentum equations are considered in the inviscid two-fluid model. Surface 
tension is also neglected since it only acts on small scales, while the waves determining the flow 
structure in pipe flows are usually of long wavelength. The gas phase is assumed 
incompressible, as the Mach number of the gas phase is usually very low for stratified flow. 
Hence, the governing equations are as follows: 


|;(«/)+^( w /«/) = 0 0 ) 

at ox 
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p g dx 


— g cos 9H g 


da, . 

-^--a g gsmO 
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where t and x are the respective time and axial coordinates, a is the volume fraction, u is the 
velocity, p is the pressure, p is the density, H is the hydraulic depth, g is the gravitational 
acceleration, (3 is the angle of inclination of the pipe; the subscripts / and g denote the liquid and 
gas, respectively, and the subscript /' denotes the interface. 

11.2. Computational Procedure 

The governing Eqs. (1-4) are solved iteratively. The basic procedure is to solve the continuity 
equation of liquid, Eq. (1) for liquid volume fraction, and liquid and gas phase momentum 
equations are used to obtain liquid and gas phase velocity. However without an appropriate 
pressure fluid, the new calculated liquid volume fraction and liquid and gas velocity may not 
satisfy the gas continuity equation. The key lies in the handling of pressure term, which cannot 
be explicitly calculated from these four governing equations. However, a governing equation for 
pressure can be easily constructed by combining Eq. (1) and Eq. (2) to a total mass 
constraining equation 


£(“,«>^0w)=O (5) 


Substituting the liquid and gas momentum equations into (5) yields a constraint on pressure. 
SIMPLE type of pressure correction scheme [10] is then used. 


A finite volume method is employed to discretize the governing equations. A staggered grid 
(Fig. 2) is adopted to obtain a compact stencil for pressure [15]. On the staggered grids, the flow 
properties such as volume fraction, density and pressure are located at the center of the main 
control volume, and the liquid and gas velocities are located at the cell face of the main control 
volume. 


Velocity control 
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Fig. 2. Staggered grid arrangement in two-fluid model. 
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The Euler backward scheme is employed for the transient term. The discretized liquid continuity 
equation becomes 


Y t ((“/ )/> “ ( a i )p ) + ( a i U l l ~ ( a i U l )w = 0 ( 6 ) 

where the superscript 0 represents the values of last time step. The symbol P refer to the center 
of main control volume, symbols e and w refer to the east face and west face of main control 
volume respectively. The liquid velocity at the cell face is known, and the volume fraction at the 
cell face can be evaluated using various interpolation schemes. Among them, central difference 
(CDS), 1 st order upwind (FOU), 2 nd order upwind (SOU) and QUICK schemes are commonly 
used. Eq. (2) for the gas phase is similarly discretized. 

The liquid momentum equation is integrated on the velocity control volume. Using similar 
notations, one obtains 

“7 ) P - (a ,11, )° p )+ (u, \ (a,u, \ - (u, ) w ( a,u , ) w = 

(a ) (7) 

— ^ (p w -Pe)+ ((«/ ),, - («/ )e) H iS cos 0 - A x(a, ) p g sin 0 

Pi 


It is important to note that the interpolation schemes used in Eq. (7) must be exactly the same 
as those in Eq. (6) in order to reduce the dissipation and dispersion errors. 

The gas phase momentum equation is similarly treated, 


Ax 

At 


(k M * l - ( a s u 8 X )+ k i ( a s u s i - k l ( a s u g l = 

(a ) 

(Pw -Pe)+ ((«/ )w - («/ )e ) H gS C0S 6 ~ ), 


gsin# 


( 8 ) 


For the pressure correction scheme, Eq. (5) is integrated across the main control volume. The 
discretized equation is 


) e - ( a g U g ) w + ( a i U l ) e ~ ( a i U l )w = 0 ( 9 ) 

Because Eq. (5) is obtained by combining Eq. (1) and Eq. (2), the discretization of Eq. (9) 
should be exactly the same as that of Eqs. (2, 6). The final pressure equation is obtained by 
substituting two momentum equations, Eqs. (7, 8) to Eq. (9). 

11.3. Characteristics and Ill-Posedness 

Eqs. (1-4) form a system of 1st order PDEs and characteristic roots, X, of the system can be 
found. If X’s are real, the system is hyperbolic. Complex roots imply an elliptic system which 
causes the two-fluid model system to become illposed because only initial conditions can be 
specified in the temporal direction. Any infinitesimal disturbance will cause the waves to grow 
exponentially without bound. 

The characteristic roots of Eqs. (1-4) are 
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When g = 0 , Eq. (10) can have real roots only if A = u v = u , . Otherwise, the two-fluid model is 
ill-posed (Gidaspow, 1974). If g ^ 0 , the real roots (or well-posedness) requirement gives 
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Pl Pg g sin 6 


g J 


a 


(11) 


Eq. (11) gives the critical value for the two phase slip velocity beyond which the system 
becomes ill-posed. The two-fluid model well-posedness condition is exactly the same as from 
the IKH analysis on the two-fluid model [7], 

11.4. Von Neumann Stability Analysis for Various Discretization Schemes 

Von Neumann stability analysis is commonly used for analyzing the stability of a finite difference 
scheme. 


Ui-l Ui-o.5 Ui+o.5 

I t ) > 

X 

Pi-1 Pi Pi+0.5 

Otj_l Cti Oti+o.5 

Fig. 3. Grids index number in staggered grid for von Neumann stability analysis. 

In this derivation, the 1 st order upwind (FOU) scheme is used as an example, and both liquid 
and gas velocities are assumed positive. Discretization of Eq. (6) using FOU leads to 

At *'^ ^ + (k k " ( l ‘ l )<4 ( a i )m)=° ( 12 ) 

Splitting the variables into base value and disturbances, the linearized equation for the 
disturbance a, is 


Mi _Sp\ ( rj/ y^ _ fa y^ )_ Ui ((^ y _ fa y ^ | = ( 


where “ A ” denotes disturbance values. Expressing disturbances as 

(a,)" =sE n e Ikx 


(13) 


(14) 
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( 15 ) 


(i,):=s,E"e- 
(«,); =s,E'e 


(16) 


where E is a common amplitude factor, k is the wavenumber, and / represents imaginary unit. 
Eq. (13) is simplified to 


Ax 

At 


(l - G 1 )+ u, (l - e 1 * j] + e,a, (e^ - e~^)= 0 


where G is the amplification factor: 


G = 


E n 

Tjn-l 


(17) 


(18) 


and ^is phase angle: 


(j) = k • Ax 


(19) 


defined over [0,;r] which represents all the resolvable wave components in the computational 
domain for the given grid. Short waves correspond to the region near ^ = n 

The wave growth equation for the gas phase mass conservation equation is similarly obtained: 

f Az 




I ~ 2 ^ I n 

-s g a g \e 2 -e 2 )=0 
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For the liquid momentum equation, Eq. (10) is discretized with the FOU scheme to 
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Pi 
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( 21 ) 


which is subsequently linearized and simplified with the aid of the liquid mass conservation 
equation, 


Ax 

At 


(pl ((«/ )i + I - («Z )”4 1+ AM/ ((«/ )” + i - («/ )/-I ) 
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( 22 ) 


For the gas phase, the velocity disturbance is governed by 
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( 23 ) 
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The pressure term can be canceled by combining Eqs. (22-23), 
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Eqs. (17, 20, 24) can be written in the form of an amplification matrix. Non-trivial solutions for 
( e,e ,e^f exist only when the determinant of the matrix is zero. Hence, 


a{fj~ x ) 2 + b{(j ~ l )+ c = 0 


(25) 
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where CFL are Courant number 


CFL, =—u. and CFL = — w 

; Ax ; g Ax g 


(27) 


The values of A(^) are given in Table 1. From Eq. (25), the amplification factor can be easily 
found, 


G = 


2a 


-bLylb 2 -4 ~ac 


(28) 


Stability requires |G| < 1 for all </> . 
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Table. 1. A (ft) for different numerical schemes. 


Scheme 

A(«0 

1 st order upwind 

X-e 1 * 

Central difference 

2 

2 nd order upwind 

3 - 4e~^ +e- 21 * 

2 

QUICK 

le 1 * +3-7e^ +e~ 21 * 
8 


11.5. Inviscid Kelvin-Helmholtz (IKH) Analysis 

IKH analysis provides a stability condition for the two-fluid model as well as useful information 
on the growth rate of disturbance in the two-fluid model. Eqs. (1-4) are linearized and 
substituted for the perturbed liquid volume fraction, liquid and gas phase velocities, and 
pressure in the form of e exp(/(<«r - Ax)) in which e is the amplitude, ryis the angular frequency, 
and k is the wavenumber. The following system is obtained for the disturbance amplitude: 
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The dispersion relation is obtained as: 
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(30) 


where c is the wave speed. It is note that the negative imaginary part of co determines growth 
rate of perturbation wave in two-fluid model. Eq. (30) is identical to Eq. (10), only with . being 
replaced by c. Details of IKH analysis can be found in [7], 
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III. Results and Discussions 


111.1. Computational Stability Assessment based on von Neumann Stability Analysis 

Comparison of stability of the 1 st order upwind, 2 nd order upwind, central difference, and QUICK 
schemes are conducted first for flow conditions before near, and after the instability. 

It is well known that for ordinary convection-diffusion equations the 1st order upwind scheme is 
less accurate with high numerical diffusion, and high order schemes, such as SOU, CDS, and 
QUICK, have lower numerical diffusion. In this study, for illustration purposes, water and air are 
considered and the pipe diameter is taken to be 0.078m. The computational domain is 1m long, 
the grid number is N = 200. The pipe incline angle is 6 is 0. The unperturbed values are: 
a , = 0.5 , u, = 1 mis, u g = 17 ml 5 and CFL t = 0.1 . Stability condition based on Eq. (10) or (30) 

for the above parameters is AU = u g - u, < 16.0768 . Thus, the two-fluid model for this condition 

is well-posed. It serves as an ideal testing case to assess the performances of various schemes 
since the system is quite close to being ill-posed. There are two values of G given by Eq. (28) 
and the larger one determines the instability; so that only the larger growth rate is used here. 

Fig. 4 compares the amplification factors G of four discretization schemes. The solid line is the 
amplification factor by IKH analysis (G = 1). The dotted line is for the CDS scheme. It is slightly 
lower than one but quite close to one with a slight diffusion error at high wavenumber range. 
This implies the CDS is an ideal scheme to compute the two-fluid model. The dashed line is for 
the FOU scheme which possesses excessive numerical diffusion at high k. Furthermore, G > 1 
at low k. Thus the computation using FOU is unstable under this flow condition. The dashed and 
dotted line is for the SOU scheme. Although SOU is regarded as a better scheme than FOU 
with less numerical diffusion, its performance in the two-fluid model is very poor. For large k, the 
numerical diffusion of SOU is much larger than that of FOU. For small k, the amplification factor 
of SOU is much larger than that of FOU. Dashed double dotted line is the amplification factor of 
the QUICK scheme. Its numerical diffusion at high k is lower than that of FOU and SOU, but it is 
still much larger than that of CDS. At small k, G is slightly larger than 1 indicating that QUICK is 
unstable as well. The reason that the amplification factor of CDS is close to the analytical 
amplification factor is probably due to a lack of 2nd order diffusion error and low dispersion error. 
Compared with SOU, numerical diffusion of FOU scheme is much higher and dispersion is 
slightly lower. Overall performance of FOU is better than that of SOU which suggests that the 
dispersion error in the two-fluid model is much more important than it is in the simple 
convection-diffusion equation. The interpolation of QUICK is essentially linear interpolation with 
the upwind correction. Therefore its numerical diffusion and stability is worse than that of CDS, 
but better than that of FOU and SOU. 

Next, the effect of the slip velocity AU = u g - u u on the numerical stability is discussed. Fig. 5 
shows the amplification factors of the CDS scheme for a range of values of AU . 
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Fig. 6. Amplification factor of FOU scheme at different A U =u g -u l . N = 200 , 6 = 0 , a, = 0.5 , 

u , = \m/ s , and CFL l =0.1. 

When A/7 is smaller than that given by the IKH stability, the amplification factors of all the 
harmonics in the computational domain are less than one. However, if A U is higher than the 
IKH stability criteria, the two-fluid model is analytically ill-posed, and G for small k exceeds one, 
as shown by the curve for A U = 16.1 mis, in Fig. 5. From numerical results, a neutral stability 
condition of CDS is found to be near A U = 16.0773m/ s m for the condition used in Fig. 5, which 
is quite close to IKH stability condition of 16.0768 m/s. As A U further increases, G increases 
too. The range of unstable harmonic wavenumber becomes wide. The amplification factor of 
CDS scheme matches that of IKH only at very low wavenumber. In the high k range, numerical 
damping causes G to be much lower than one. 

Fig. 6 shows the amplification factors of the FOU scheme at different values of A/7 . Unlike the 
CDS scheme, there is no significant change of G when At/ varies. Numerical results indicate 
that the neutral stability for the condition shown in Fig. 6 is A/7 = 14.772 mis, which is much 
lower than the analytical value of 16.0768/m . The behavior of SOU and QUICK scheme is 
similar to the FOU. The stability condition for SOU is A/7 = 13.73 mis and for QUICK it is 
A/7 = 16.03m/ s . 

Fig. 7 shows the effect of the liquid velocity on the amplification factors of CDS with 
AU = \6m/s and A//Ax = 0.1 . For u,=0.0\m/s and u,=0.lm/s , G decreases 

monotonically with the phase angle. Damping appears at high k. When u, increases, G at high 
k range rises significantly, leaving a high damping saddle at the intermediate k range. On the 
other hand, if A/7, is constant, CFL g /CFL, is much larger than one when u, is small so that it 
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is hard to keep both CFL, and CFL g in the moderate range, which is essential to the 
computational stability. 


Fig. 8 shows the effect of u, on G for the FOU scheme with AU = 1 6m I s and At/ Ax = 0.1. The 
behavior of FOU is much different with that of CDS. When I u is small, most harmonics are 
unstable. For a larger u , , excessive numerical diffusion makes the computations stable. 
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Fig. 8. Amplification factor of FOU scheme at different u r N = 200, A U = 16 mis 

6 = 0 = 0.5 , and At/ Ax = 0.1. 
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Fig. 9. Amplification factor of CDS scheme at different At/ Ax . N = 200 , 
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Fig. 10. Amplification factor of FOU scheme at different At/ Ax . N = 200 , 

Uj = 1 mis, AU = 1 6m Is, 0 = 0, and a, = 0.5 . 


Figs. 9-10 show the effect of At /Ax on G for the CDS and FOU schemes. Both show increasing 
numerical diffusion with increasing At/ Ax resulting in a decrease in G. This can be explained 
by examining Eq. (26c), where the last term is related to the gravitational effect. It is well known 
that gravity stabilizes the flow. Thus increasing At/ Ax computationally helps the stability. 

IV. Conclusions 

Numerical instability for the incompressible two-fluid model near the ill-posed condition is 
investigated for various discretization schemes, while the pressure correction method is used to 
obtain the pressure. The von Neumann stability analysis is carried out to obtain the amplification 
factor of a small disturbance in the discretized system. The central difference scheme has the 
best stability characteristics in handling the two-fluid model, followed by the QUICK scheme. It is 
quite interesting to note that the excessive numerical diffusion in the 1 st order upwind scheme 
seems to promote the numerical instability in comparison with the central difference scheme. 
Despite its nominal 2 nd order accuracy and popularity, the 2 nd order upwind scheme is much 
more unstable than the 1 st order upwind scheme for solving two-fluid model equations. Different 
discretization schemes for the convection term with varying degrees of the numerical diffusion 
and dispersion cannot cause a delay in the stability; they often promote instability in the two- 
fluid model. 

The analytically predicted wave amplitude growth is also compared with that obtained from 
carefully implemented computations using various discretization schemes for the convection 
term. Excellent agreement between the numerical results and the predicted results is obtained 
for the growth of the wave amplitude and the dominant wavenumber when the computation 
becomes unstable. The relation between the computational instability and the ill-posedness is 
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discussed. In the presence of the smallamplitude long-wave disturbance, whose amplitude is 
much larger than the machine round-off error, the growth of the disturbance exactly matches the 
prediction of the von Neumann stability analysis when the computational stability condition is 
violated. In the meantime, a shorter wave emerges from the machine round-off error, and 
eventually dominates the entire disturbance, which causes the computation to blow up. This 
computational instability is widely interpreted as the result of ill-posedness of the two-fluid model. 
The results of the present study suggest that the computational instability is largely the property 
of the discretized two-fluid model and is strongly affected by the inherent ill-posedness of the 
two-fluid model differential equations. Introduction of numerical diffusion and/or dispersion can 
significantly change the instability of the discretized system; however, such steps often yield 
unfavorable computational results. For solving two-fluid models, central difference is 
recommended since it is much more accurate and dependable than other schemes investigated. 

V. Patents, Publications, Presentations and Students from Research 

Publications 

1. Liao, J., Mei, R., Klausner, J.F. 2005, “A Study on Numerical Instability of Inviscid Two- 
Fluid Model Near Ill-Posedness Condition”, Proceedings of the ASME Summer Heat 
Transfer Conference, San Francisco, CA. 

Students from Research 

Jun Liao: Recived PhD. Employed with Westinghouse Electric Company, Monroeville, PA. 

VI. Acknowledgements 

This Work was supported by NASA Glenn Research Center under contract (NAG3-2750) 
and NASA Kennedy Space Center. 


VII. References 

1. Wallis, G.B., 1969, One-Dimensional Two-Phase Flow, McGraw-Hill, New York. 

2. Ishii, M., 1975, Thermo-Fluid Dynamic Theory of Two Phase Flow, Eyrolles, Paris. 

3. Gidaspow, D., 1974, “Modeling of Two Phase Flow”, Proceedings of the Fifth 
International Heat Transfer Conference, VII, pp. 163. 

4. Jones, A.V., and Prosperettii, A., 1985, “On the Stability of First-Order Differential 
Models for Two-Phase Flow Prediction”, International Journal of Multiphase Flow, 11, 
pp. 133-148. 

5. Song, J.H., and Ishii, M., 2000, “The Well-Posedness of Incompressible One- 
dimensional Two-Fluid Model”, International Journal of Heat and Mass Transfer, 43, 
pp. 2221-2231. 

6. Issa, R.I., and Kempf, M.H.W., 2003, “Simulation of Slug Flow in Horizontal and Nearly 
Horizontal Pipes with the Two-Fluid Model”, International Journal of Multiphase Flow, 29, 
pp. 69-95. 

7. Barnea, D., and Taitel, Y., 1994, “Interfacial and Structural Stability of Separated Flow”, 
International Journal of Multiphase Flow, 20, Suppl., pp. 387-414. 

8. Brauner, N., and Maron, D.M., 1992, “Stability Analysis of Stratified Liquid-Liquid Flow”, 
International Journal of Multiphase Flow, 18, pp. 103-121. 

9. Chan, A.M.C., and Banerjee, S., 1981, “Refilling and Rewetting of a Hot Horizontal Tube 
part II: Structure of a Two-Fluid Model”, Journal of Heat Transfer, 103, pp. 287- 292. 

10. Patankar, S.V., 1980, “Numerical Heat Transfer and Fluid Flow,” McGraw-Hill, New 
York. 


NASA/CR— 2008-21 5440/PART3 


138 



Il.lssa, R.I., and Woodburn, P.J., 1998, “Numerical Prediction of Instabilities and Slug 
Formation in Horizontal Two-Phase Flows”, 3rd International Conference on Multiphase 
Flow, ICMF98, Lyon, France. 

12. Lyczkowski, R.W., Gidaspow, D., Solbrig, C.W., and Hughes, E.D., 1978, 
“Characteristics and Stability Analyses of Transient One-Dimensional Two-Phase Flow 
Equations and Their Finite Difference Approximations”, Nuclear Science and 
Engineering, 66, pp. 378-396. 

13. Stewart, B.H., 1979, Stability of Two-Phase Flow Calculation Using Two-Fluid Models, 
Journal of Computational Physics, 33, pp. 259-270. 

14. Ohkawa, T., and Tomiyama, A., 1995, “Applicability of High-Order Upwind Difference 
Methods to the Two-Fluid Model”, Advances in Multiphase Flow, Elsevier Science, 
pp. 227-240. 

15. Ferziger, J.H., and Peric, M., 1996, Computational Methods for Fluid Dynamics, 
Springer, Berlin. 

16. Hirsch, C., 1988, “Numerical Computation of Internal and External Flows, volume I: 
Fundamentals of Numerical Discretization”, John Wiley & Sons, New York. 

17. Liao, J., 2005, “Modeling Two-Phase Transport During Cryogenic Chilldown in a 
Pipeline”, Ph.D Dissertation, University of Florida. 


NASA/CR— 2008-21 5440/PART3 


139 



Sub-Task 3: Numerical Investigation of Cryogenic Fluid Transport in Horizontal Pipelines 
During Turbulent Chilldown Process 

Abstract 

A numerical model has been proposed to predict the chilldown history of pipe wall temperature 
during the transport of cryogenic fluid in horizontal pipeline under unsteady high mass flow rate 
regime. The model then is compared with existing experimental data for verification. There is 
very good agreement between simulation results and experiment data of Jackson et al. (2005) 
measured at the four locations on the pipelines. The proposed model has captured important 
heat transfer features between phases and between liquid and pipe wall, and also provides a 
very inexpensive tool for predicting thermal history during chilldown process. 

Nomenclature 

c p heat capacity 

D inner pipe diameter 

d thickness of pipe wall 

g gravity 

h heat transfer coefficient 

h poo i pool boiling heat transfer coefficient 

h conv convection boiling heat transfer coefficient 

h fg latent heat 

k thermal conductivity 

p pressure 

Nu Nusselt number 

Re Reynolds number 

S suppression factor 

T temperature 

T Lei Leidenfrost temperature 

T crit transition temperature from nucleate boiling to convection heat transfer 

T 0 ambient temperature 

t time 

U velocity 

Greek Symbols 

z, r, and cp cylindrical coordinates 
a liquid volume fraction 

5 vapor film thickness 

p viscosity 

6 dimensionless temperature 

a liquid surface tension 

Subscripts 

i and o inner and outer pipe 
/ liquid 

v vapor 

w wall 

sat saturated 

fb film boiling 

nb nucleate boiling 
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I. Introduction 


1.1. Literature review of experimental and numerical studies on cryogenic Chilldown 

Cryogenic chilldown process is encountered in many engineering applications. One of the most 
important things is the transport of cryogenic propellants in internal fuel tanks of space shuttle in 
space shuttle launch facilities through a complex pipeline system. To maintain liquid fuel 
entering the fuel tanks, a cryogenic chilldown prior to the filling is required to cool the pipe wall 
so that its temperature goes down to the saturated temperature of the core cryogenic liquid fuel. 

Cryogenic chilldown process has intrinsic complication due to the interactions among two 
phases of the cryogenic fluid and the solid. It is a complicated problem due to the coupling 
between fluid dynamics and heat transfer. The studies of cryogenic chilldown process emerged 
during 60’s decade parallel to the development of space technology. The so-called “space race” 
between developed countries during these times has promoted researches involving cryogenic 
engineering, both experimentally and computationally. 

Several experimental studies of cryogenic fluids were conducted in the early of 1960s. Burke et 
al. [1] studied chill down in various stainless-steel lines by flowing liquid nitrogen. A model to 
predict chill down time was developed by treating the entire line as a single control volume. 
They have assumed an infinite heat transfer rate from the cryogenic fluid to the pipe wall. The 
effects of flow regimes on the heat transfer rate were neglected. This system has provided a 
simple estimate of chill down time but lacked accuracy due to the averaging of fluid properties 
and flow rates over the chill down time. The experimental results reported by Bronson et al. [2] 
shown that stratified flow is prevalent during cryogenic chilldown system. 

Flow regimes and heat transfer regimes in horizontal pipe during cryogenic chilldown were also 
studied by Chi and Vetere [3]. Several flow regimes were identified: single vapor phase, mist 
flow, slug flow, annular flow, bubbly flow, single liquid phase. Heat transfer regimes were also 
identified as single phase vapor convection, film boiling, nucleate boiling, and single phase 
liquid convection. Graham et al. [4] correlated heat transfer coefficient and pressure drop with 
Martinelli number based on their experimental data. 

Chi [5] developed a one-dimensional model for energy equations of the liquid and the wall, 
based on the film boiling heat transfer between the wall and fluid. This model of the chilldown 
was developed under the assumptions of constant flow rate, constant heat transfer coefficient, 
constant fluid properties, homogeneous fluid flow, and film-boiling-dominated heat transfer. An 
empirical equation for predicting chilldown time and temperature was proposed. The assumption 
of constant flow rate is not very realistic because usually the transfer line inlet and exit 
pressures are set while the flow rate may vary greatly according to the flow condition. The 
assumptions of constant heat transfer coefficient and constant fluid properties are highly 
restrictive, but can provide useful estimates of temperature and chill down time. 

Steward [6] numerically developed a homogeneous flow model for cryogenic chilldown using 
finite difference formulation. The model treated the cryogenic fluid as a homogeneous mixture. 
The continuity, momentum and energy equations were solved simultaneously to obtained 
density, pressure and temperature of the mixture. Various heat transfer regimes were 
considered: film boiling, nucleate boiling and single-phase convection. Heat transfer coefficients 
were determined using superposition of single-phase forced convection correlations and pool 
boiling correlations for both nucleate and film boiling. Correct modeling of flow film boiling is 
especially crucial with cryogenic hydrogen because film boiling occurs to as low as a tube wall 
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to fluid saturation temperature difference of approximately 1 1 K resulting in a large portion of the 
chill down occurring under film boiling. 

The stratified flow regime, which is the prevalent flow regime in horizontal chilldown, was first 
studied by Chan and Banerjee [7-9]. They developed a separated flow model for the simulation 
of the quenching by a stratified flow in a hot horizontal pipe. Both phases were modeled using 
one-dimensional mass and momentum conservation equations. The wall temperature was 
computed using a two-dimensional transient heat conduction equation. Their prediction for the 
wall temperature agreed well with their experimental results. Although significant progress was 
made in handling the momentum equation, the heat transfer correlations employed were not as 
advanced. 

Recently, Velat et al. [10] systematically conducted experiments involving cryogenic chilldown 
with liquid nitrogen in horizontal pipe. Their studies included a visual recording the chilldown 
process in a transparent Pyrex pipe (used to identify flow regimes and heat transfer regimes), 
collecting temperature histories at different positions of the wall in chilldown, and recording 
pressure drop along the pipe. Chung et al. [11] conducted similar study with nitrogen chilldown 
at relatively low mass flux and provided data needed to assess various heat transfer coefficients. 

Liao et al. [12] developed an analytical model for film boiling heat transfer for stratified laminar 
flow of cryogenic fluid. By analytically solving the mass, momentum and energy equations in the 
vapor film layer, an analytic solution of the vapor film thickness is obtained. Then, the film 
boiling heat transfer in the vapor film layer is obtained simply by the ratio of the thermal 
conductivity and the thickness of the vapor film. This model is then tested with experimental 
data of Chung et al. [11] and some agreement is achieved. 

Jackson et al. [13] conducted experiment with liquid nitrogen to measure the heat transfer 
coefficients in nucleate boiling regime during chilldown process. The recorded temperature 
history of the pipe wall and numerical study of the unsteady heat transfer through pipe wall are 
used to determine heat transfer coefficients which then are compared with well-known nucleate 
flow boiling heat transfer correlations of Gungor and Winterton, Kandlikar, Muller-Steinhagen, 
etc. However, due to the complication of the unsteady, high mass flux cryogenic flow, the 
predicted heat transfer coefficients obtained by this inverse technique do not agree well with 
current correlations. 

1.2. Heat transfer regimes of cryogenic chilldown 

A typical chilldown process involving several heat transfer regimes is shown in Fig. 1. Due to the 
large difference between temperature of the pipe wall (typically at laboratory temperature) and 
the cryogenic temperature of the entering liquid, film boiling occurs in which the liquid phase is 
supported on a thin vapor layer adjacent to the wall. As shown in Fig. 1, near the liquid front is 
the thin vapor film boiling. The knowledge of the heat transfer in the film boiling is limited 
because it is not the major interest of industrial applications and more importantly the 
experiment carrying out to capture the behavior at this stage is extremely difficult due to the high 
difference in temperature as mentioned above. 
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Fig. 1 . Heat transfer regimes of cryogenic chilldown process inside a horizontal pipe. 

There are several works from Bromley [14], Dougall and Rohsenow [15] and Laverty and 
Rohsenow [16] that showed experimental results on the vapor film boiling on vertical surfaces. 
Film boiling on horizontal surfaces was first studied also by Bromley [14] and the Bromley’s 
correlation was widely used. Then, Breen and Westwater [17] modified Bromley’s correlation to 
take into account the various changes in diameter of the horizontal tubes. Empirical correlations 
for cryogenic film boiling were also proposed by Hendrick et al. [18-19], Ellerbrock et al. [20], 
von Glahn [21]. Giarratano and Smith [22] gave a more detailed assessment of these 
correlations. All the correlations proposed above were constructed based on experiment for 
steady state cryogenic film boiling. Therefore, their validity for transient chilldown applications is 
still suspicious. 

When the pipe wall chills down further, film boiling ceases and nucleate boiling occurs in which 
the temperature of the interior pipe wall falls below the Leidenfrost temperature. In this stage, 
the liquid phase comes into contact with the pipe wall. For simplification, it is often assumed that 
the switch from film boiling to the nucleate boiling is immediate, not passing through a transition 
boiling regime. The position from film boiling transitioning to the nucleate boiling is often called 
the rewetting front as the liquid phase stars to contact with the pipe wall. Normally, we utilize the 
Leidenfrost temperature as an indicator of the transition from the film boiling to the nucleate 
boiling. 

There are numerous studies on the forced convection boiling heat transfer. A general correlation 
for saturated boiling was first introduced by Chen [23], In this work, he developed a two-phase 
flow boiling correlation where the heat transfer coefficient was divided into a nucleate pool 
boiling coefficient and a bulk convective coefficient. Bennett and Chen [24] modified the Chen’s 
correlation using a larger data set. They implemented a new correlation for the nucleate pool 
boiling suppression factor. The two-phase enhancement factor was still a function of the 
Lockhart-Martinelli two-phase multiplier. Gungor and Winterton [25] modified Chen’s correlation 
and made the extension to subcooled boiling. They also introduced the enhancement and 
suppression factors for macro-convective heat transfer. Their expression for the pool boiling 
convective heat transfer coefficient included reduced pressure, molecular weight, and heat flux. 
Kutateladze [26] and Steiner [27] also provided correlations for cryogenic fluids in pool boiling 
and forced convection boiling. These correlations are applicable for cryogenic fluids since they 
are developed directly based on cryogenic conditions. 
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As the wall temperature decreases further over time or the wall superheat drops to a certain 
range, nucleate boiling is suppressed and heat transfer regime is governed by single phase 
forced convection [28], [32]. 

II. Formulation 


To gain the fundamental insight into the thermal interaction between the cryogenic fluid and 
solid pipe wall and to predict chilldown process varies in time, various heat transfer coefficients 
must be evaluated. For heat transfer in solid pipe wall, we need to solve the diffusion equation 
to obtain temperature field and then obtain heat flux from inner pipe wall into the cryogenic liquid. 
The heat transfer between the cryogenic fluid and the wall is modeled using different heat 
transfer correlations depending on the operating heat transfer regime at a given location. These 
issues will be addressed in details in the following sections. 

11.1 . Heat transfer in solid pipe wall 


The unsteady three-dimensional heat conduction equation in cylindrical coordinates in the solid 
pipe wall is as follows [29]. 


Ps C 


p,s 


dT 

dt 
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r dr 


f 




r d(p\r dtp ) 


( 1 ) 


where p s is the material density of the solid pipe wall, c ps is the specific heat capacity and 
k s is the thermal conductivity of the solid pipe wall, respectively. T is the temperature, t is the 
time variable and ( r,<p,z ) denote the cylindrical coordinates. Since the first term on the RHS of 
Eq. (1) is small compared with the other terms [29], it is therefore neglected. Eq. (1) becomes: 
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Equation (2) can be further non-dimensionalized by using the following non-dimensional 
parameters: 
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where T sat is the saturated temperature of the cryogenic fluid, T w is the temperature on the 
outer wall of the pipe, d is the thickness of the pipe wall. The subscript “0” denotes the 
properties at the reference temperature T 0 . Then, Eq. (2) is transferred to: 
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The boundary condition for the outer pipe wall is written as follows. 
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(4a) 




where h out is the heat transfer coefficient between the outer pipe wall with the surroundings. 

It is usual that h out is very small as one wants to reach the insulated boundary condition. h out is 

obtained by experimental measurement. T amb is the ambient temperature and T out is the 

temperature at the outer surface of the pipe. 

The boundary condition for the inner pipe wall is as follows. 

K-^ = h Jt (T„-T,.) (4b) 


where h Jh is the heat transfer coefficient between the inner pipe wall with the vapor film due 
to high wall superheat during the cryogenic chilldown process. The subscript “fb” stands for “film 
boiling”. T in is the temperature at the inner surface of the pipe. h fb is a very important quantity 

since film boiling plays a major role in heat transfer between cryogenic liquid and solid pipe wall. 
Detailed computational approach for estimating h jb , which is the major task of this work, is 

given in details in section 11.2.2.1. 

11.2. Heat transfer between cryogenic fluid and solid pipe wall 

The cryogenic fluid entering the pipe firstly is in liquid phase. However, due to the high 
superheat from pipe wall, in which the pipe temperature usually is in room temperature, part of 
the cryogenic liquid vaporizes into vapor phase. Hence, the cryogenic fluid inside the pipe after 
initial setup comprises both liquid and vapor phases. Therefore, we need to calculate both heat 
transfer regimes between vapor phase and solid pipe wall and between liquid phase and solid 
pipe wall. 

The heat transfer coefficient between vapor phase and solid pipe wall is computed via available 
correlation. However for the heat transfer between cryogenic liquid and solid pipe wall, as we 
already addressed above, comprises three stages: film boiling, nucleate boiling and forced 
convection boiling. Numerical investigation of heat transfer coefficient for nucleate boiling and 
forced convection boiling is not the main interest in this work; hence, we use available 
correlations with some modifications to compute these heat transfer regimes. The main issue 
here is to analytically and numerically investigate the heat transfer coefficient during film boiling 
regime since currently there exists no suitable correlation for the present problem. 

11.2. 1 Heat transfer between cryogenic vapor and solid pipe wall 

Heat transfer between cryogenic vapor and solid pipe wall can be predicted by assuming that 
the vapor flow is fully developed convection vapor flow, i.e., we neglect the cryogenic liquid 
droplets entrained in the vapor flow which is due to the fact that the vapor flow is moving faster 
than the liquid flow. The heat transfer coefficient for single phase turbulent forced convection 
vapor flow is adopted from Dittus-Boelter’s correlation [30] and is a given as follows. 
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Re v and Pr v is the vapor Reynolds number and Prandtl number respectively. 


Re. = 


»vA> 

M v 


( 6 ) 


Pr, = 


c u 

p,vr*v 


(7) 


Here u v is the vapor velocity which is measured from the experiment, p v is the vapor viscosity, 
k v is the vapor thermal conductivity, c is the heat capacity of vapor phase and D v is the 
hydraulic diameter of the vapor phase. 


11.2.2. Heat transfer between cryogenic liquid and solid pipe wall 

As we mentioned above, heat transfer regimes between cryogenic liquid and solid wall consists 
of film boiling, nucleate boiling and forced convective boiling. The transition from one to another 
depends on many parameters and flow conditions. For simplicity, we exploit the Leidenfrost 
temperature as the transition point to determine the boiling heat transfer regime. If wall 
temperature is higher than the Leidenfrost temperature, film boiling dominates. If the wall 
temperature is between Leidenfrost temperature and a critical temperature T crit , nucleate boiling 

occurs. If the wall temperature is below T crit , single phase forced convective heat transfer is 

assumed. It is hard to tell exactly what Leidenfrost temperature and critical temperature T crit for 

different types of cryogenic fluid because they highly depend on various flow conditions. In this 
work, they are extracted from experimental results where we observe the sudden changes in 
measured solid wall temperature over time. 

11.2.2. 1. Film boiling heat transfer 

a. Introduction 


Due to the high wall superheat in cryogenic chilldown process as the cryogenic liquid enters the 
pipes in room temperature, film boiling heat transfer regime play a very important role in terms 
of both time span and fraction amount of heat removed from the solid pipe wall. It is to be noted 
again that there is no specific correlation available for this film boiling regime that has an 
intrinsic high wall superheat. If one wants to use existing conventional correlations, one needs 
to make significant modifications because those correlations were established on non-cryogenic 
conditions. It must also be noticed that correlations that were obtained for steady flow film 
boiling are not appropriate for the chilldown process due to the inherent unsteadiness of the 
present problem. 

There are several correlations for film boiling heat transfer in the literature that were based on 
the analyses of the thin vapor film boundary layer and the stability of the thin vapor film. Bromley 
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[14] and Breen and Westwater [17] introduced correlations for film boiling on the outer surface 
of a hot tube. Ferderking and Clark [31] and Carey [32] introduced correlations for film boiling on 
the surface of a sphere. 

Recently, Liao et al. [12] developed an analytical model for film boiling heat transfer of stratified 
flow during cryogenic chilldown process. This model was obtained by analytically solving the 
mass, momentum and energy equations in the vapor film layer. Since this model is only applied 
for stratified laminar cryogenic flow and with further assumptions, all the conservative equations 
can be solved analytically and an analytic solution of the vapor film thickness is obtained. Then, 
the film boiling heat transfer in the vapor film layer is obtain simply by the ratio of the thermal 
conductivity and the thickness of the vapor film due to the linear temperature profile. This model 
is the first attempt for solving cryogenic film boiling. However, it is only applicable for low mass 
flux flow regime. For flow with high mass flux, this model is no longer valid and we need a 
suitable model to take into account the effect of high flow rate which is more prevalent in 
industrial processes. 

With high mass flow rate, the liquid and vapor film flows are supposed to be in turbulent regime. 
Then, the problem turns out to be convective heat transfer in turbulent boundary layer problem 
and solving for this type of problem is the main issue of this work. 


i 



Fig. 2. Schematic of film boiling model. 

The schematic plot of the model for film boiling inside a horizontal pipe is shown in Fig. 2. The 
bulk liquid is near the bottom of the pipe. Between the liquid and the solid pipe wall is a thin 
vapor film as a result of the occurrence of film boiling. Due to buoyancy effect, the vapor in the 
film flows upward along the azimuthal direction. However, due to the high wall superheat, the 
vapor film is produced to separate the cryogenic liquid from touching the solid pipe wall. The 
balance between these two physical phenomena leads to presence of the thin vapor film and 
heat is transfer from the solid wall to the liquid through this thin vapor film. 
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There are certain assumptions that can be made to simplify the analysis; however, the physical 
insight of the problem is still remained. These assumptions are as follows. 

- The liquid flow is incompressible and turbulent. 

- The vapor film thickness 5 is small compared to the solid pipe radius/? and is a function of 
position only: 8 = 8{cp) . 

- Phase change occurs at the liquid-vapor interface. 


b. Formulation 

The continuity, momentum and energy equations for that thin vapor film boundary layer are as 
follows. 


du dv 
— + — = 0 

dx dy 


du du 1 dp d 2 u 

u h v — = + v v — -- gsmcp 

dx dy p v dx dy 2 


( 8 ) 
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The momentum in 
hydrostatic pressure 


dp 

/-direction yields — = 0 , then p = p(x ) = const along y-direction. The 

dy 

can be evaluated by: 

P = P 0 +P,gR(cos(p-cos(p 0 ) (11) 


Then: 


dP 

— = ~Pig smp 
dx 


Using the following well-known Reynolds decompositions 


( 12 ) 


u = U + u' 
v = V + v' 
p = P + p' 

Then the governing equations become as follows: 

dU dV_ 
dx dy 


(13) 
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In Eq. (15) p v — is the molecular shear stress acting in x-direction, - pu'v' = ps M — is the 

dy dy 

turbulent stress, and s M is the eddy viscosity. 

dT dT 

In Eq. (16) a v — is the heat flux due to molecular motion, - pc v'T' = pc e H — is the 

dy dy 

turbulent heat flux, and s H is eddy thermal diffusivity. 


Because the turbulent diffusion across the thin vapor film is typically much greater than 
convection in azimuthal direction, the convective terms in momentum and energy equations are 
neglected. Thus the further simplified momentum and energy equations and the corresponding 
boundary conditions are given as follows. 

Simplified momentum equation: 


d 

dy 


(A + £ m ) 


dU^ 
dy A 


= (A ~P,)g^(p 


Boundary conditions: 

1 . At wall: no-slip condition: 

2. At the liquid /vapor interface: no frictional force: 
Simplified energy equation: 


y = 0,U = 0 

a dU n 
v = o, = 0 

dy 
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dy 
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(A + £ h ) 
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cry 

dy J 


= 0 


Boundary conditions: 

1 . At wall: y = 0,T = T W 

2. At the liquid /vapor interface: y = S,T = T sat 


(17) 


(18a) 

(18b) 


(19) 


(20a) 

(20b) 


At this stage, in order to solve the simplified momentum and energy equations, we need a 
model to approximate the eddy viscosity and eddy thermal diffusivity. The well-known Kays’s 
model [33] is employed. 
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c. Analytical solution 
Velocity field 

From Eq. (17), let: 


Mx) = ( Pi ~ P v )g sin (p = ( A - p v )g sin 


\Rj 


Then, integrate Eq. (17): 


dU - A(x) C, 

— = — — y + 1 

dy jU v +s M (y ) ju v +s M (y ) 

At_y = 0 : y + =0 then from Eq. (21): s M {y = 0) = 0 , hence from Eq. (26): 


dU _C L 

dy v=0 Pv 


Then from Eq. (22): 
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At y = 8 : exploiting Eq. (18b), Eq. (26) yields: 


dU Q -A(x) Si C, 
dy y=s Mv + £ m (d) /* v + e M (S) 
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And then: 


C, = A( x)S (30) 

With this expression for integration constant C x , we can obtain the following expression for 

u in Eq. (22) and in Eq. (26) as follows. 

dy 


* 

u 


lA(x)8 

V Pv 


(31) 


dU 

dy 


A(x)8 1 

V v +e M (y) L 




= /O0 


(32) 


At this stage, we can integrate Eq. (32) with information from Eq. (21), (22), (25), and (31) to 
obtain the velocity field. 

Temperature field 

Integrate Eq. (19) we obtain: 


dT _ C 2 

dy a v +s H (y) ( 33 ) 

We then integrate Eq. (33) from 0 to 5 and using the two boundary conditions given by 
Eq. (20a) and (20b), we obtain the expression for C 2 as follows. 


T -T 
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r dy 
a v +e H (y) 


Then, plug Eq. (34) into (33) yields: 
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(35) 


At this stage, we can integrate Eq. (35) with information from Eq. (21), (22), (23), (24), (25), (31), 
and (34) to obtain the temperature field. 

Constraint on vapor film thickness 

The energy and mass balance on the vapor film requires the following condition: 
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or: 
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Where U is the mean velocity and given as follows by definition: 


u =\ i u “y 


(37) 


And from Eq. (35): 


dT_ 

dy 


C, 


y=8 


CC V + £ H (<5) 


(38) 


Eq. (37) is the key for this analysis since it is the constraint for the vapor film thickness. Once $ 
is obtained, the film boiling heat transfer coefficient is computed as follows. 




q_ 
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, dT 
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(39) 


Where A T sat =T W - T sat is the wall superheat. 

Interestingly, we may recover the analytical result of Liao et al. [12] for the film boiling heat 
transfer coefficient. By plugging Eq. (38) into Eq. (39) and using Eq. (34), the resulting equation 
is given as follows. 


hjb ~ 


dy 


a v +e H {y) 
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(«v £ H (8)) 


For Liao et al. model, there is no eddy viscosity as well as eddy thermal diffusivity since they 
consider the laminar flow regime (or low mass flux flow regime), hence by letting 

e H (y) = e H (y) = 0 in Ec l- ( 4 °) y' elds: 
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Since a v is constant and 



= 8 , we recover Liao et al. result [12]. 
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( 42 ) 



K 
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By recovering the analytical result from previous work of Liao et al. [12], we have just shown 
that our analysis is not only consistent but it also is the general case of the low mass flux flow 
regime. 

11.2.2.2. Nucleate boiling heat transfer 

In the nucleate boiling stage, there are several correlations are available in the literature and are 
commonly used, such as Chen’s [23], Gungor and Winterton’s [25], Kutateladze’s [26]. Liao et al. 
[12] suggested that a modified Kutateladze’s correlation is suitable for the nucleate boiling heat 
transfer of cryogenic liquid. 

In this work, a modified Kutateladze’s correlation [26] with a suppression factor is used. 

h nb= h conv+ Sxh pool ( 43 ) 


Here h conv is given by the Dittus-Boelter’s correlation [30] for turbulent fully developed pipe flow 
and h pool is the pool nucleate boiling heat transfer coefficient. They are given as follows. 
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Re, and Pr ; is the vapor Reynolds number and Prandtl number respectively. 


Re, = 


u,D, 

Pi 
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c _pjP, 
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(46) 

(47) 


Here u, is the liquid velocity which is measured from the experiment, //, is the liquid viscosity, 
k, is the liquid thermal conductivity, c , is the heat capacity of liquid phase, D, is the hydraulic 

diameter of the vapor phase, p is the pressure, <r is the surface tension and A T mt is the wall 
superheat. 

11.2.2.3. Forced convective heat transfer 
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When the wall superheat drops to a certain range, all the nucleate sites are suppressed [32], 
then single phase forced convection heat transfer dominates. The heat transfer coefficient in 
this case can be predicted by using Dittus-Boelter’s correlation [30] for turbulent flow regime and 
given in Eq. (44). 

III. Results and Discussions 

111.1 . Experiment of Jackson et al. 

In experiment by Jackson et al. [13], liquid nitrogen was stored in a 1580 kPa dewar and used 
as cryogen. The dewar pressure was the driving force for the liquid nitrogen flow. Upon leaving 
the dewar, nitrogen was cooled by passing through a heat exchanger prior to entering the 
vacuum jacketed concentric Pyrex visual test section. The test section has inner diameter of 
12.57 mm and the thickness is 1.65 mm. The flow structure was captured by a high resolution 
CCD camera system. 

The flow then entered the heat transfer measurement section located at 70 cm downstream of 
the visual test section. In this section there were two thermocouples measure the temperature of 
the fluid at the inlet and outlet. The inner and outer pipe wall temperature in several 
circumferential positions were measured by thermocouples and recorded on a computer. The 
thermocouple arrangements and dimensions of the test section are shown in Fig. 3. 

Experiments were carried out at room temperature and atmosphere pressure. The mass fluxes 
ranged from 66 to 500 kg/m 2 s, vapor qualities varied from 0.004 to 1 and ambient heat fluxes 
ranged from 716 W/m 2 to 995 W/m 2 . 

The Leidenfrost temperature for the nitrogen in this experiment was observed around 130K; 
hence the temperature when the film boiling ends and nucleate boiling starts is set at 130K. The 
transition temperature at which the nucleate boiling switches to single phase convection heat 
transfer was observed around 105K based on experimental results. 



Fig. 3. Thermocouple positions in the circumference of the Pyrex pipe. 

III. 2. Computational Procedure 

Finite Volume Method is employed in computing the thermal field in the solid pipe wall. The 
computational domain consists of 40 grids in the radial direction and 40 grids in the azimuthal 


NASA/CR— 2008-21 5440/PART3 


154 



direction. The grid resolution is sufficient in capturing the thermal behavior. The computational 
procedure is as follows. 

First, we need to read in the experimental data including the following information: 

• Dimension of test section 

• Solid and fluid properties at initial stage 

• Raw vapor velocity data 

• Raw liquid velocity data 

• Raw liquid height data 

• Raw liquid pressure data 

Then, for each time step, we need to implement the computational procedure as shown in Fig. 4. 
It can be interpreted as follows. 
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Fig. 4. Flow chart of computational procedure for each time step. 
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• First we guess a temperature profile at the inner solid pipe wall. 

• Then, by comparing the inner temperature at every location in the azimuthal direction 
with the Leidenfrost and critical temperature, we obtain the corresponding heat 
transfer coefficient for each type of boiling regimes. 

• For film boiling heat transfer, we need to iterate to obtain the correct vapor film 
thickness in order to get film boiling heat transfer coefficient. First we guess a value 
of vapor film thickness, then following the procedure addressed in section 2.2.1, we 
obtain the velocity and temperature field throughout the vapor film layer. The iteration 
for vapor film thickness stops when the mass and energy balance satisfies. 

• After getting all the heat transfer coefficients at the inner pipe wall, it is obviously that 
we have all the information needed for boundary condition at the inner pipe. We also 
have the information of ambient temperature as well as the experimental heat 
transfer coefficient of the surroundings; we can compute the boundary condition at 
the outer pipe wall. 

• We then solve for the temperature field in the solid pipe wall by finite volume method 
with the above two boundary conditions. The system of equation is solved efficiently 
by Thomas algorithm. 

• The iteration for inner wall temperature stops when the difference between the new 
obtained inner wall temperature with the guessed one less than a chosen criterion. 

• Record the temperature information and then we move to the next time step. 

III. 3. Comparison of pipe wall temperature 

Since the major effort is in the modeling of film boiling stage, only numerical results in this 
boiling regime are shown. Following is some comparison between the numerical results 
obtained by exploiting the above developed model with the experimental results. 

The mass flow rate for the first experiment ranges from 80 to 100 kg/m 2 s. The numerical 
simulation started at 20s since as the experiment needs that amount of time for starting-up. 

For the top pipe wall, the computed temperature is less than that measured from experiment 
about 10 to 20% as shown in Fig. 5. 

For the side pipe wall, the discrepancy between computed temperature and experiment result is 
still large but less than that for the top pipe wall as shown in Fig. 6. 

For the bottom pipe wall, we can observe a very good agreement between the numerical result 
and experimental result for the temperature history as shown in Fig. 7. 

The mass flow rate for the second experiment ranges from 130 to 200 kg/m 2 s. The numerical 
simulation started at 20s since as the experiment needs that amount of time for starting-up. The 
same observation for top, side and bottom walls can be concluded from Figs. 8, 9, and 10. 

From above observations, we see that for the top and the side pipe wall, the predicted 
temperature is not so good. It can be explained by some of the following reasons. First reason 
may be due to the error in recording experiment results. Remember that the input for our 
computation is the raw data from experiment, hence, the numerical results depends a lot on 
those data. The second reason may be the correlations we are using are not quite suitable or 
we need to modify them in order to take into account the effect of high mass flow rate conditions. 
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However, the predicted temperature at the bottom pipe wall agrees quite well with the 
experimental results. It shows that our model for the film boiling heat transfer is a very good one 
and is also reliable. Since this is the high mass flux regime, we expect the cryogenic flows to be 
in turbulent regime and that we build the model to take into account the turbulent effects which 
manifest via eddy viscosity and eddy thermal diffusivity. 

Besides, we are also able to compute the thin vapor film thickness that varies for each location 
in the azimuthal direction for each time step. That information is very useful since it is 
impossible for us to measure and record the vapor film thickness by conducting experiment, 
especially under such unsteady and high mass flow rate regime. 


Experiment vs Numerical Temperature - Top Pipe Wall 



Fig. 5. Comparison between experimental and numerical results for temperature at top wall. 

Mass-flux ranging from 80-100 kg/m 2 s. 
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Experiment vs Numerical Temperature - Side Pipe Wall 



Fig. 6. Comparison between experimental and numerical results for temperature at side wall. 

Mass-flux ranging from 80-100 kg/m 2 s. 


Experiment vs Numerical Temperature - Bottom Pipe Wall 



Fig. 7. Comparison between experimental and numerical results for temperature at bottom wall. 

Mass-flux ranging from 80-100 kg/m 2 s. 
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Experiment vs Numerical Temperature - Top Pipe Wall 



Fig. 8. Comparison between experimental and numerical results for temperature at top wall. 

Mass-flux ranging from 130-200 kg/m 2 s. 


Experiment vs Numerical Temperature - Left Side Wall 



Fig. 9. Comparison between experimental and numerical results for temperature at side wall. 

Mass-flux ranging from 130-200 kg/m 2 s. 
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Experiment vs Numerical Temperature - Bottom Pipe Wall 



Fig. 10. Comparison between experimental and numerical results for temperature at bottom wall. 

Mass-flux ranging from 130-200 kg/m 2 s. 


IV. Conclusions 

A computational model has been developed for understanding the heat transfer 
mechanisms during cryogenic chilldown process. With this model, one can predict the thermal 
interactions among cryogenic liquid phases and between phases and solid pipe wall. Especially, 
a model for simulation of film boiling heat transfer is developed based on transport equations in 
the vapor film layer and the effect of high mass flow rate. This model is tested with experimental 
results and the agreement between the two is quite acceptable. 
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Abstract 

For many industrial, medical and space technologies, cryogenic fluids play indispensable roles. 
An integral part of the cryogenic transport processes is the chilldown of the system components 
during initial applications. In this paper, we report experimental results for a chilldown process 
that is involved with the unsteady two-phase vapor-liquid flow and boiling heat transfer of the 
cryogen coupled with the transient heat conduction inside pipe walls. We have provided 
fundamental understanding on the physics of the two-phase flow and boiling heat transfer 
during cryogenic quenching through experimental observation, measurement and analysis. 
Based on the temperature measurement of the tube wall, the terrestrial cryogenic chilldown 
process is divided into three stages of film boiling, nucleate boiling and single-phase convection 
that bears a close similarity to the conventional pool boiling process. In earth gravity, cooling 
rate is non-uniform circumferentially due to a stratified flow pattern that gives rise to more 
cooling on the bottom wall by liquid filaments. In microgravity, there is no stratified flow and the 
absence of the gravitational force sends liquid filaments to the central core and replaces them 
by low thermal conductivity vapor that significantly reduces the heat transfer from the wall. Thus, 
the chilldown process is axisymmetric, but longer in microgravity. 

Keywords: Cryogenics: Boiling heat transfer; Two-phase flow; Chilldown; Microgravity 

1. Introduction 

Cryogenic fluids are widely used in industrial, aerospace, cryosurgery systems and so on. For 
example, liquid hydrogen is used in industrial applications such as metal processing, plate glass 
production, fat and oil hardening, semiconductor manufacturing, and pharmaceutical and 
chemical manufacturing. In these systems, proper transport, handling and storage of cryogenic 
fluids are of great importance. However, the chilldown or quenching process which initiates 
cryogenic fluids transport is complicated, involving unsteady two-phase heat and mass transfer, 
and has not been fully understood. 

1.1 Role of Cryogenics in Space Exploration 

The extension of human space exploration from a low earth orbit to a high earth orbit, then to 
Moon, Mars, and possibly asteroids and moons of other planets is one of NASA’s biggest 
challenges for this new millennium. Integral to this is the effective, affordable, and reliable 
supply of cryogenic fluids. The efficient and safe utilization of cryogenic fluids in thermal 
management, power and propulsion, and life support systems of a spacecraft during space 
missions involves the transport, handling, and storage of these fluids in terrestrial and reduced 
gravities. 

Chilldown process is the inevitable initial stage during cryogenic fluid transport. Due to their low 
boiling points, boiling and two-phase flows are encountered in most of the cryogenic operations. 
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The complexity of the problem results from the intricate interaction of the fluid dynamics and 
heat transfer, especially when phase-change (boiling and condensation) is involved. Because of 
the large stratification in densities between the liquid and vapor phases, the reduced gravity 
condition in space would strongly change the terrestrial flow patterns and accordingly affect the 
momentum and energy transport characteristics. Therefore, boiling and two-phase flow behave 
quite differently when the gravity levels are varied. The uncertainties about the flow pattern, 
pressure drop and heat transfer characteristics pose a severe design concern. For example, 
there is considerable disagreement over the chilldown heat fluxes and whether a unique 
rewetting temperature exists [1-2]. Different definitions of rewetting temperatures are reported in 
chilldown researches: Leidenfrost temperature, minimum film boiling temperature, quenching 
temperature from thermocouple observations and the temperature at which critical heat flux 
occurs [1]. For similar experimental observations, quiet different explanations were also 
suggested by different researchers. For example, it was reported that the rewetting velocity 
increased with increasing inlet flow rate, given the same initial wall temperature [3-5], Duffey 
and Porthouse [4] suggested that the flow rate effect resulted from increasing the wet side heat 
transfer coefficient with higher inlet flow. This improves the rate of axial heat conduction and 
hence leads to a faster rewetting rate. Thompson [5], however, argued that the inlet flow rate 
affects precooling on the dry side rather than the heat transfer in the wet side. Additionally, 
because of the experimental difficulties, in general, there is very little heat transfer data for 
cryogenic flow boiling in reduced gravity. A representative work is reported by Adham- 
Khodaparast et al. [6]. 

The proposed research will focus on addressing specific fundamental and engineering issues 
related to the microgravity two-phase flow and boiling heat transfer of cryogenic fluids that 
require well-designed and meaningful experimentation. In this paper, liquid nitrogen chilldown 
process is investigated and divided into several stages based on the flow image and 
temperature data. The outcome of the research will provide fundamental understanding on the 
transport physics of cryogenic boiling and two-phase flows in a reduced gravity. 

2. Background and Literature Review 

2.1 Boiling Curve 

A boiling curve shows the relationship between the heat flux that the heater supplies to the 
boiling fluid and the heater surface temperature. According to the typical boiling curve as shown 
in Fig. 1, a cryogenic chilldown (quenching) process usually starts from point E and then goes 
towards point D in the film boiling regime as the wall temperature decreases. Point D is called 
the Leidenfrost point which signifies the minimum heater temperature required for the film 
boiling. For the film boiling process, the wall is so hot that liquid will vaporize before reaching the 
heater surface which causes the heater to be always in contact with vapor. When cooling 
beyond the Leidenfrost point, if a constant heat flux heater were used, then the boiling would 
shift from film to nucleate boiling (somewhere between points A and B) directly with a 
substantial decrease in the wall temperature because the transition boiling is an unstable 
process. If the heater wall temperature can be controlled independently, then the boiling 
process will proceed from D to C in the transition boiling mode. 
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Free convection Nucleats Transition Film 



Fig. 1. Typical boiling curve. 


2.2 The Effects of Gravity 

Because of the differences in density and inertia, the two phases in the flow are usually non- 
uniformly distributed across the pipe under the terrestrial condition. The absence of gravity has 
important effects on flow patterns, pressure drop and heat transfer. Surface-tension-induced 
forces and surface phenomena are likely to be much more important in outer space than they 
are on earth. Actually, all flow-pattern-specific phenomena will be influenced by the gravity level 
[7], As an example, in one of our own experiments, for the same mass flow rates, Fig. 2 
compares the flow patterns under both terrestrial and microgravity conditions and the difference 
is significant as the liquid filament is on the bottom and suspended in the middle of the tube for 
terrestrial and microgravity conditions, respectively. 
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Fig. 2. Gravity effect on flow patterns. A) Flow pattern in 1-g test. B) Flow pattern in 
microgravity test. 


Momentum and energy transport in two-phase flow and boiling under microgravity is a relatively 
new field. This is a complex problem for the scientific community to solve. The rewards are rich, 
however, as new ways to propel future space missions and more intensive space activities have 
been planned. Currently, there exists a limited quantity of literature, data, and test facilities to 
even begin to answer these questions. Microgravity cryogenic boiling research was initiated 
more than fifty years ago, but the progress has been limited. The following summarizes the 
accomplishments. 

2.3 Terrestrial Cryogenic Boiling and Two-Phase Flow 

Numerous studies of cryogenic boiling in the one-g environment were conducted in the 1950’s 
and 60’s. Brentari et al.’s work [8] is a comprehensive review of the experimental studies and 
heat transfer correlations. For oxygen, nitrogen, hydrogen and helium, it was found that for pool 
boiling, the Kutateladze [9] correlation had the greatest reliability for nucleate boiling while the 
Breen and Westwater [10] correlation was best for film boiling. Maximum nucleate flux data 
were reasonably well predicted by the Kutateladze correlation. Although these correlations were 
selected as the best available, neither has particularly good agreement with experimental data. 
For the case of forced convection boiling, Brentari et al [8] reported that no correlation was 
found to be distinctly better. Some simple predictive methods were found to work as well as 
more complex schemes. In all boiling cases, it was questioned as to whether or not the 
predictive correlations include all of the significant variables that influence the boiling process. In 
particular, it was suggested that more detailed and better controlled experiments are needed 
and that more attention to surface and geometry effects is required. 

Another comprehensive review of cryogenic boiling heat transfer addressing hydrogen, nitrogen 
and oxygen is given by Seader et al. [11], It was reported that nucleate pool boiling results 
cannot be correlated by a single line but cover a range of temperature difference for a given 
heat flux. The spread is attributed to surface condition and geometry, and orientation. Maximum 
heat flux can be reduced by about 50% when going from normal-g to near zero-g. Seader et al 
[11] reported a fair amount of data for film pool boiling. Film boiling heat flux is reduced 
considerably at near zero-g conditions. Only a very limited amount of data is available for 
subcooled or saturated forced convection boiling and few conclusions were drawn. Relatively 
recent correlations have been published for normal-g saturated flow boiling of cryogens using 
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the Convection, Boiling and Froude numbers as correlating parameters [12-15], For example, 
Van Dresar et al. [15] experimentally studied the near-horizontal two-phase flow of nitrogen and 
hydrogen. Unlike most of the other works which based on turbulent liquid flow, their work 
focused on laminar liquid flow conditions and the results for low mass and heat flux flow were 
correlated with Froude number. 

2.4 Reduced Gravity Cryogenic Boiling and Two-Phase Flow 

In general, there is little heat transfer data for cryogenic flow boiling in reduced gravity. We were 
able to find just two reports. Adham-Khodaparast et al. [6] have investigated the flow film boiling 
during quenching of R-113 on a hot flat surface. They reported lower heat transfer rates during 
microgravity as compared to normal gravity and contributed that to thickening of the vapor layer. 
The wall superheat and the surface heat flux at the onset of rewetting and maximum heat flux 
were found to increase with the inlet liquid subcooling, mass flux and gravity level. The effect of 
gravity was determined to be more important for low flow rates and less relevant for high flow 
rates. Antar and Collins [16] reported flow visualization and measurements for flow film boiling 
of liquid nitrogen in tubes on board KC-135 aircraft. They were particularly interested in the 
vapor/liquid flow pattern and the thermal characteristics. They identified a new vapor/liquid flow 
pattern that is unique in low gravity. They also observed that a sputtering leading core followed 
by a liquid filament annular flow pattern. This new flow pattern is composed of a long and 
connected liquid column that is flowing in the center of the tube and is surround by a thick vapor 
layer. The vapor annulus that separates the liquid filament from the wall is much thicker than 
that observed in the terrestrial experiment. They attributed the filamentary flow to the lack of 
difference in the speed of vapor and liquid phases. On the heat transfer side, they reported that 
the quench process is delayed in low gravity and the tube wall cooling rate was diminished 
under microgravity conditions. 

3. Experimental System 

To investigate the cryogenic chilldown process, a liquid nitrogen two-phase flow experimental 
facility has been designed, fabricated and tested under both terrestrial and microgravity 
conditions. A drop tower is used to provide the microgravity condition. 

The experimental system is designed as a once-through flow pass using either gravity or motor- 
driven bellows as flow generator. Fig. 3 shows the schematic of the experimental system, which 
locates in two side-by-side aluminum cubicles and is fabricated for both terrestrial and 
microgravity experiments. The experimental system mainly consists of a nitrogen tank, a motor- 
driven bellows, test section inlet portion, test section, test section outlet portion, vacuum jacket, 
vacuum pump, data acquisition system, lighting and video system. A photographic view of the 
apparatus is shown in Fig. 4. Each component in the system is described separately below. 

3.1 Flow Delivery Mechanisms 

We have used two different flow delivery systems, one is a gravity-driven flow and the other is a 
constant flow rate bellows-driven flow. In the gravity-driven system, an insulated reservoir is 
used to generate the flow. By reviewing the recorded flow images, the mass flux in the gravity 
flow is estimated to be between 18-23 kg/m 2 s which provides a liquid entering flow velocity 
around 3-4 cm/s. 
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Standard Rig 




Fig. 4. Experimental apparatus located in cubicles. 
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Fig. 5. Cryogenic flow driven system. 


Table 1. Bellows-driven flow delivery system parameters 


Motor speed 

Liquid velocity entering the 

Mass flux 

Approximated test 

(rpm) 

test section (cm/s) 

kg/m 2 s 

duration 
limit (minute) 

5 

0.446 

3.606 

15 

10 

0.891 

7.205 

7.5 

15 

1.337 

10.811 

5 


In the bellow-driven system, we followed the traditional method to control cryogenic flow, where 
the nitrogen flow is generated by a motor-driven stainless steel bellows (Fig. 5). The basic idea 
is using a constant speed motor to pull a moving plate which is attached to the bottom of a 
bellows filled with liquid nitrogen. The bellows is submerged in liquid nitrogen inside a nitrogen 
tank. During the experiment, the state of the flow at the exit of the bellows is saturated liquid. As 
the bellows being compressed at a constant speed, a constant volumetric flow rate can be 
achieved, and the corresponding mass flow rate can be calculated. Details of the flow delivery 
component are given in Yuan et al. [17], Table 1 shows the three motor speeds used in the 
experiment and the corresponding mass flux, flow velocity and maximum time of operation. It is 
noted that the bellows-driven flows are considered lower flow rates while the gravity driven flows 
are higher flow rates. 

3. 2 Test Section 

The test section is a Pyrex glass tube of 25.4 cm long. The ID and OD of the test section are 
11.1 mm and 15.8 mm, respectively. The test section inlet and outlet are stainless steel tubes. 
At both ends of the test section, stainless steel adaptors and Teflon ferrules are used to connect 
the test section to the test section inlet and outlet portion. 
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There are nine drilled holes of approximately 2mm depth in the test section. The diameter of 
each hole is 1 mm. A total of 16 type-T thermocouples are placed on the test section, nine are 
embedded very close to the inner surface through drilled holes at three downstream cross- 
sections. At each cross-section, three thermocouples are located circumferentially at equal 
separation distance. The other seven thermocouples are used to measure the outside wall 
temperatures at two cross-sections, also located circumferentially at equal separation distance. 
The test section can be rotated along its axis before being fastened at two ends. Fig. 6 sketches 
the test section and the thermocouple locations in one of the tests. Fig. 7 is a photograph of the 
test tube with attached thermocouples. 




V \ 


Section 4 10,12,13 Section 2 7,8,9 


J 


33**in 

Fig. 6. Test section and thermocouple locations. 



Fig. 7. Photograph of cryogenic two-phase flow test apparatus. 

The test section inlet, test section and test section outlet are enclosed in a vacuum jacket built 
from stainless steel vacuum components. Two transparent quartz windows in the vacuum jacket 
enable the observation and record of the two-phase flow regimes inside the test section. The 
diameter of each window is 7.62 cm. A ceramic sealed vacuum feed-through flange is used to 
connect the thermocouple wires from the vacuum side to the air side. The vacuum is maintained 
by a potable vacuum pump during the experiments. A CCD camera (CV-730 from Motion 
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Analysis Inc.) faces one of the quartz windows to record flow images, while lighting is provided 
by a fluorescent light at the other widow. 

3.3 Drop Tower Apparatus 

The microgravity environment was provided by the University of Florida drop tower. The drop 
tower is 5-story high and is equipped with an airbag for deceleration. Fig. 8 shows the overall 
drop tower schematic. At the current drop height of 15.25 m, UF’s drop tower produces 
1.7 sec of free fall. The experiment is actually located in a drag shield which prevents the 
experiment from the aerodynamic drag. As a result, the microgravity level is measured between 
10' 5 to 1CT 4 g. Details of the drop tower system and its operation are given in Yuan et al. [17], 


Winch 



Fig. 8. A schematic of the drop tower system. 

3.4 Data Acquisition Component and Experimental Uncertainties 

A data acquisition system is built for recording temperatures and flow images during 
experiments. Type-T thermocouples (Omega) are used for temperature measurement. The 
thermocouples are wired to a screw terminal board and then connected to a 16-channel 
thermocouple board (PCI-DAS-TC from Measurement Computing) plugged into the PCI slot of a 
computer. All the thermocouples are tested and calibrated with boiling nitrogen prior to the 
chilldown experiments. A Labview program is developed to read the temperature measurements 
to the computer. Video images are monitored and recorded by connecting the CCD camera to a 
frame grabber board. A commercial software records the flow images and also shows the real- 
time images on the computer screen. 
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For the bellows-driven system, the uncertainty of mass flux mainly comes from accuracy of the 
driven-motor speed, and it is calculated that the relative error for mass flux measurement is 
6.88%. The type-T thermocouples used for temperature measurement have the uncertainty of 
+0.5 °C declared by the manufacturer. For highly transient process like chilldown the response 
time of the thermocouples is also important. To get quick response, the tip style of the 
thermocouples is chosen as exposed and the wire diameter of the thermocouples is used as 
smaller as possible. With the wire diameter of 0.25 mm the responding time is less than 
0.2 second according to the chart given by the manufacturer. 

Another uncertainty source of temperature measurement comes from the data acquisition (DAQ) 
system. The DAQ board for temperature measurement has programmable gain ranges and A/D 
pacing, and accepts all the thermocouple types. The accuracy of the measurement depends on 
the gain, the sample rate and the thermocouple type. The uncertainty of type-T thermocouple is 
±0.9 °C for worst case from the product specification. For current experiment, the gain is set at 
400 and the sample rate is about 60 Hz. It is found that the uncertainty for current settings is 
about +0.3 °C. The relative error for mass flux measurement is +6.88%. 

4. Experimental Results and Discussion 

During the “pipe chilldown” experiment, cryogenic transport pipes experience a fast cooling on 
the walls due to an initially large temperature difference (liquid nitrogen at 77K and wall at 300K). 
The purpose of the chilldown experiment is to provide an overall characterization of the system 
that includes both flow visualization and heat transfer measurements, and also to offer flow 
patters and heat transfer characteristics during the pipe chilldown process that takes place in 
many cryogenic applications in a room temperature environment. The following discussions are 
classified into two categories: low flow rates (bellows-driven) and high flow rates (gravity driven). 
In the beginning of the chilldown experiment, liquid nitrogen enters the Pyrex glass tube that is 
at the room temperature (25 °C). The bulk velocities of the entering liquid flows are estimated in 
the range of 0.5 to 4 cm/s with both gravity-driven and bellow-driven techniques. With this flow 
range, the role of gravity is measured by the dimensionless Froude number, Fr. which defined 
as the ratio of flow inertia to gravitational force, U 2 /gD, where U is the entering bulk velocity, g is 
the gravitational acceleration constant and D is the tube diameter. If the Fr number is much 
larger than unity, then the flow inertia is dominant that renders the gravity insensitive. On the 
other hand, if the Fr number is much less than unity, then gravity is dominant. The range of the 
Fr number corresponding to the flow velocities between 0.5 and 4 cm/s is from 1.83x10" 4 to 
1.47x1 O' 2 . Since the Fr number is less than unity by two orders of magnitude, it is safe to 
conclude that gravity plays an important role for both gravity-driven and bellow-driven flows, but 
the gravity is totally dominant for the bellows-driven case. Because the glass tube is unheated, 
the tube wall is going through a transient cooling process. In the following, the measured wall 
temperature history is presented first and then flow visualization pictures are given to provide 
the correlations between the wall temperature and flow pattern. 

4.1 Gravity-Driven Flows - High Flow Rate Case 

4.1.1 Boiling Phenomenon during Quenching 
A. Terrestrial Gravity Case 

For the chilldown process, the temperature histories of the tube walls were measured by the 
thermocouples at three downstream locations (near inlet, middle and near outlet). The results 
are given in Fig. 9. Based on the characteristics of the temperature history curves, the chilldown 
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process is generally divided into three stages: (i) Film boiling, (ii) Nucleate boiling and (iii) 
Single-phase convection. 

(i) Film boiling. At the beginning of the process, the wall temperature of the test section is 
near the room temperature, which is much higher than the temperature of the liquid nitrogen 
and above the Leidenfrost pint, therefore film boiling is taking place on the tube wall surface. 
Based on the flow visualization presented in the next section, the initial quenching stage is 
associated with a hot wall which causes the two-phase flow to take the form of a vapor core with 
scattered small liquid chunks (Dispersed Flow Film Boiling, DFFB). In general, the temperatures 
measured by thermocouples that are located in the bottom of the tube (Channels 12, 15, 9, 6, 
and 3) are lower than those from thermocouples embedded in the upper portion. This is mainly 
due to the gravity effects that bring liquid fragments cooling to the bottom walls. 





Time (s) 

c 


Fig. 9. Temperature profiles at three downstream locations during the chilldown process. A) 
Temperature profile of section 1 (near inlet), B) Temperature profile of section 2 (middle point), 
C) Temperature profile of section 3 (near outlet). 
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(ii) Nucleate boiling. A distinctive character of the bottom wall temperature history is the 
sharp increase in the slope at the wall temperature of about 190 K. This phenomenon is seen at 
all three downstream locations. The physical explanation is based on the transition of boiling 
regimes from film boiling to transition boiling and then to nucleate boiling as mentioned earlier 
and shown in Fig. 1 . Even though, the transition boiling regime is very short. 

(iii) Single-phase convection. Finally when the wall is chilled down enough such that the 
wall temperature (<140 K, point A in Fig. 1) is too low to support nucleate boiling, the heat 
transfer switches to single-phase forced convective boiling 

(iv) Chilldown boiling curve. A 2-D transient conduction model is used to calculate heat 
transfer data from the transient temperature profiles [17]. In this model, energy balance is 
performed locally on a control volume of the tube wall at inside wall thermocouple location. The 
change in stored heat in the control volume is equated to the heat transported to the fluid and 
heat transferred by conduction, minus losses to the environment. The error of this model is 
evaluated to be in the range of 5-1 1 %. 

Based on the temperature data, the heat transfer coefficient history for the bottom wall at 
section 1 is estimated using the temperature data and plotted in Fig. 10. Left side of Fig. 10 
shows the calculated bottom wall heat flux as a function of bottom wall temperature at the outlet 
cross-section, while the right side represents the corresponding temperature profiles re-plotted 
from Fig. 9C. 


Wall Temperature (K) 



Fig. 10. Estimated heat transfer coeficient. 

The shape of the heat flux on the bottom of the tube is similar to the boiling curve from steady- 

state pool boiling experiments. This suggests that the chilldown process may share many 

common features with the pool boiling process. Following the method to characterize different 
heat transfer mechanisms in pool boiling experiments, a local maximum or critical heat flux 
(CHF) and a local minimum heat flux q" mm are used to divide the chilldown heat transfer 

into three stages, which are film boiling, transition boiling and nucleate boiling, as shown in 

Fig. 10 similar to that defined in Fig. 1 for pool boiling. 
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Fig. 1 1 . Bottom wall heat flux at the outlet cross section as a function of time. 

Initially, the wall temperature is very high and above the Leidenfrost point, so the liquid phase 
can not contact the tube wall (non-wettable). As a result, liquid nitrogen evaporates intensively 
when entering the test section; a vapor film blanket will cover the wall surface and separate the 
liquid from contacting the wall, the two-phase flow is therefore in the film boiling state. As the 
wall temperature decreases during quenching and drops below the Leidenfrost temperature, the 
liquid begins to contact the wall; the heat transfer is in the transition boiling regime, which is 
characterized by an increasing wall heat flux with a decreasing wall superheat that is contrary to 
the trends in the film boiling region and nucleate boiling region. After passing through the CHF 
point the heat transfer mechanism then enters the nucleate boiling regime. The calculated 
bottom wall heat flux as a function of time is shown in Fig. 11. It is obvious that the time of 
transition boiling is very short as compared with the other two boiling stages. 

The similarity between the chilldown boiling curve and the pool boiling curve naturally leads one 
to compare the chilldown data with pool boiling correlations. The comparison between the two 
turning points, namely the minimum heat flux and the CHF, is given below. For a steady state 
film boiling, the correlation developed by Zuber [18] is widely used to predict the minimum heat 
flux: 


q" = Ch, p 

“mm Iv r v 


<rg{p,~P v ) 

(a -a ) 2 


( 1 ) 


Here.cr is surface tension; gis gravitational acceleration; h lv is the latent heat of vaporization, 
while C , a constant, was suggested by various researchers as 0.177 [18], 0.13 [19], or 0.09 [20], 
The resulting g" in is then 13.0 kW/m 2 , 9.6 kW/m 2 , or 6.6 kW/m 2 , respectively. In Fig. 11, the 

q" mia for chilldown is calculated as 13.3 kW/m 2 , which is slightly larger than the steady state 

prediction with C = 0.177. Based on the similarity between the CHF condition and the column 
flooding, Kutateladze [21] derived the following relation for the pool boiling CHF: 


NASA/CR— 2008-21 5440/PART3 


177 


q" cllF =0.m h ,vp v 


( 2 ) 


°g{pt-p v )y 

2 


Zuber [18] developed the identical correlation based on the analysis of Taylor and Helmholtz 
instability. For liquid nitrogen under the atmospheric pressure, Equation (2) gives a CHF value 
of 160.7 kW/m 2 . The chilldown measurement in Fig. 1 1 is only about 27% of this value. This big 
discrepancy is believed to come from the following possible sources: (i) Equation (2) was 
developed for pool boiling over a large horizontal cylinder or a sphere, while the current case is 
boiling inside a small tube, and (ii) different experimental conditions between the chilldown and 
the pool boiling processes. In the pool boiling experiment, the heat supplied to the fluid is 
maintained by a heater, while in chilldown tests it comes from the stored heat in the tube wall. In 
the film boiling region, the heat flux is generally small, therefore the tube wall can maintain a 
near constant heat flux condition and function like the heater used in the pool boiling tests. 
However, in the transition boiling region, the stored energy in the tube wall is depleted so 
quickly that the experimental condition is very different form that of the pool boiling tests. The 
limited energy stored in the tube wall put a restriction on the value of CHF, which, therefore, is 
much less than the pool boiling data. Previous work by Bergles and Thompson [22] also 
quantitatively concluded that the differences between chilldown and steady-state boiling curves 
can be very large. 

For Eqs. (1) and (2), the minimum heat flux and the CHF, thermal properties of the wall are not 
factored in the equations. However, as seen from the above, the available heat flux to the flow is 
closely related to the energy stored in the wall. Therefore, the thermal properties, e.g., thermal 
conductivity, heat capacity, of the wall are expected to play a role in the chilldown correlations. 
For example, for a wall with higher thermal conductivity, the energy transferred to the fluid can 
be more quickly supplied by the walls, and therefore it will result in a higher value of CHF. 



Temperature (K) 

Fig. 12. Bottom wall heat fluxes at different axial locations. 

Fig. 12 shows the estimated heat fluxes in different axial locations of the tube. It is found that 
both the CHF and the minimum heat flux decrease with increasing axial distance from the inlet. 
This is also associated with an increase in the rewetting temperature. 
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In most of the previous quenching experiments, the film boiling heat flux was reported as either 
a relatively constant value along the axial direction [23;24] or that it decreases downstreams as 
the test section is chilled down [22;25]. In our experiments, however, the local heat flux first 
increases in a short period then decreases gradually. This is consistent with the transient nature 
of the experiments. Generally the increasing period is expected to be shorter at higher mass 
fluxes if the other conditions are kept the same. The mass fluxes in previous investigations were 
much larger than that in current experiments, and therefore associated with very a short time 
period of increasing heat fluxes in the film boiling regime. This might be the reason that the 
increases of the heat flux were not recorded before. 

B. Microgravity Case 

Because that the gravity-driven flow delivery system could not fit into the limited two-cubicle 
frame as shown in Fig. 4, therefore we do not have high-flow rate microgravity results. However, 
the microgravity drop tower results are provided in the low flow rate section. It is also noted that 
the gravity effects are dominant for the low flow rate case as mentioned before. But for higher 
flow rates, we can speculate that the wall temperature history curve in microgravity would be 
very similar to those of the two upper portion locations (Ch. 1 1 and 14 or 13 and 16) where the 
walls are in contact with vapor phase only. In microgravity, there would be no stratification and a 
perfect inverted annular flow would prevail, which causes the vapor to cover all the tube wall 
and leave liquid in the central core area. 

4.1.2 Two-Phase Flow Dynamics During Quenching 

A. Terrestrial Condition 

This phase of the experimentation was intended to provide not only some physical 
understanding of the two-phase flow characteristics during the chilldown process, but also 
the performance evaluation for the apparatus which includes both image and data acquisition. 



Vapor Layer 
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Liquid Layer 



Fig. 13. Two-phase flow images under 1-g condition for chilldown process. A) Initial stage 
of chilldown-vapor core with scattered liquid chunks, B) Initial stage - inverted annular 
flow, C) Transition stage, D) Transition stage, E) Cold wall - stratified flow, F) Cold wall - 
stratified flow. 


The chilldown process is divided into three stages based on the temperature of the tube wall as 
discussed above. The initial stage as shown in Fig. 13A is associated with a hot wall which 
causes the two-phase flow to take the form of a vapor core with scattered small liquid chunks. 
When the wall becomes slightly cooler, the inverted annular flow is seen as given in Fig. 13B 
where the liquid is occupying most of the tube with a vapor layer adjacent to the tube wall. For 
both cases as shown in Fig. 13A and Fig. 13B, it is believed that film boiling is the main 
mechanism (portion between D and E in Fig. 1). When the wall temperature decreases further, 
the two-phase flows are experiencing a transition from the film boiling (inverted annular flow) to 
the stratified flow. Fig. 13C and Fig. 13D show the transition flows for warmer tube walls where 
liquid chunks are starting to accumulate near the bottom side of the tube wall due to the gravity 
effects. Finally when the wall is chilled down enough, a stable stratified two phase flow was 
observed. As shown in Fig. 13E and Fig. 13F, a stable and continuous liquid layer is located on 
the bottom of the tube and a pure vapor core is flowing on top of the liquid layer. As discussed 
in the previous section that for the stratified flow here, nucleate boiling (portion between A and B 
in Fig. 1) is taking place between the liquid layer and the cold bottom wall while the film boiling 
is the mode between the vapor core and the tube wall. Therefore under the stratified flow, there 
is a circumferential temperature variation on the tube wall. 

As mentioned before, the liquid entering velocity is around 3-4 cm/s but the gravity force is still 
two orders of magnitude larger than the inertia force as mentioned previously. Therefore, the 
gravity effect overpowers the flow inertia. This is the reason why we observed the stratified flow 
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when the wall temperature drops below the Leidenfrost point that allows the wall to be wetted by 
the liquid phase. 

B. Microgravity Case 

Fig. 14 provides a series of two pictures for the chilldown experiment in a microgravity condition. 
We followed a large chunk of liquid as it traveled downstream. It is seen that basically it is an 
inverted annular flow where the vapor is filling the tube with liquid chunks in the center. The size 
of the liquid chunk is decreasing due to vaporization as it flows downstream. 



A B 

Fig. 14. Photographs of flow pattern in microgravity. 

4.2 Bellows-Driven Flows - Low Flow Rate Case 

4.2.1 Boiling Phenomenon During Quenching 

First, for lower flow rates, the gravity dominates more than the flow inertia that makes the 
cryogenic transport in microgravity to behave differently from that on earth. Also as the flow rate 
lowers, the heat transfer generally slows down that extends the film boiling period and causes 
the vapor film to thicken. Especially in microgravity, there is no possibility of a stratified flow and, 
the heat transfer is generally lower under the microgravity condition due to the replacement of 
liquid filaments by the vapor film near the bottom walls that makes the chilldown even longer. In 
this section, the heat transfer and two-phase flow in the film boiling region during chilldown is 
investigated for both terrestrial and microgravity conditions. 

A. Wall Temperature Profiles 

To measure more data points during the microgravity period, the wall temperatures are 
measured only at one cross-section near the inlet. Fig. 15 gives the typical temperature profiles 
with three different mass fluxes. 
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Fig. 16. Wall temperature response to microgravity. A) Mass flux of 3.6 kg/m 2 s. B) Mass flux of 
7.2 kg/m 2 s. C) Mass flux of 10.8 kg/m 2 s. 
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Based on Fig. 15, the general cooling rate of the tube wall is directly proportional to the mass 
flow rate of the quenching fluid. The circled temperatures in Fig. 15 approximately indicate the 
time period of microgravity condition during one drop, which includes the release of the drag 
shield, microgravity time, impact on the air bag, and deceleration period. 

To examine the detail of the wall temperature response to sudden removal of the gravity force, 
the temperature profiles are zoomed in. Figure 16 illustrates the focused results. It is observed 
that the temperature decreasing rate at the bottom wall is generally slower during the 
microgravity period, because the removal of the gravitational force will send liquid filaments to 
the central core and replace them by low thermal conductivity vapor that significantly reduces 
the heat transfer from the wall. Since the top wall transfers heat mainly by vapor convection, the 
gravity field is expected to have negligible effects on the heat transfer at the top wall. For the 
current transient experimental condition, the best way to represent the data is to calculate the 
averaged heat flux which is given in the following section. 

B. Wall Heat Flux 

The gravity effect is shown in Fig. 17, in which the ratio of bottom heat flux before drop and 
during drop is plotted for different flow rates. The heat flux at the bottom of the tube decreases 
under microgravity condition and the ratio varies from a minimum of about 0.66 to about 0.90. 
The result does not show a strong dependence on wall temperature and inlet flow rate. Two 
runs of the quenching test performed by Xu [26] reported similar ratio of 0.7 and 0.8, however, 
in Westbye et al. [23], this ratio was found to be much less and ranged from 0.15 to 0.6. 



Wall Temperature (K) 


Fig. 17. Ratio of heat flux under microgravity to 1-g condition with different flow rates and 
comparison with model prediction. 

The bottom wall is subjected to film boiling through the liquid filaments near the wall and the 
convection to the super heated vapor phase. For the film boiling part, the heat flux and the heat 

transfer coefficient are proportional to (a/ g) as suggested by Merte and Clark [27]. Assuming 

that the convection is not affected by the microgravity condition, the resulted heat flux ratio is 
shown in the dash line in Fig. 17; the gravity level in this calculation is 10' 4 g. The calculation 


NASA/CR— 2008-21 5440/PART3 


184 




results are much scattered and significantly less than the experimental values. This suggests 
that under microgravity condition, the effect of convection part may raise the heat flux, and 
therefore the total effect of the microgravity is less prominent than that in the pool boiling 
experiments. 

4.2.2 Two-Phase Flow Dynamics during Quenching 

The flow patterns before (terrestrial gravity) and during drop (microgravity) are compared in 
Fig. 18 for lower flow rates. The images on the left are taken before the drop, while the 
microgravity images are on the right. 

For all the terrestrial cases, they are all stratified flows. For Fig. 17A the mass flow rate is the 
lowest, so the liquid flow is in the form of chucks on the bottom wall. For higher flow rates in 
Figs. 17C, 17E, and 17G, the liquid flow on the bottom is in a continuous stream. Different flow 
behaviors have been recorded during the microgravity period. If the liquid phase before drop is 
dispersed liquid chunks on the bottom wall, these liquid chunks will enter the central core region 
still as chunks or stretched ligaments shown in Fig. 17B during the microgravity period. For 
higher flow rates where continuous liquid filaments are on the bottom in earth gravity, these 
filaments are lifted up and still maintain their original shapes during the microgravity period 
(Fig. 17D). In some other cases for even higher flow rates, the larger liquid filaments are broken 
and dispersed into the central region (Fig. 17F) or has a liquid-vapor core in the center and 
smaller chunks at both top and bottom (Fig. 17H). 
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Fig. 18. Two-phase flow images under both 1-g and microgravity conditions. A) 1-g case 1, B) 
Microgravity case 1, C) 1-g case 2, D) Microgravity case 2, E) 1-g case 3, F) Microgravity case 
3, G) 1-g case 4, FI) Microgravity case 4 

5. Conclusion 

This research is aimed at the fundamental understanding of the cryogenic chilldown phenomena 
inside a pipe in terrestrial gravity and microgravity. An experimental system was designed, built 
and calibrated. The microgravity environment was simulated in a 1.7-second drop tower. Flow 
patterns and wall temperature history during the chilldown were obtained for a range of different 
flow rates. The wall temperature histories were measured at three downstream locations. 
Primarily we found that the terrestrial chilldown process is divided into three stages based on 
the temperature of the tube wall. These stages are film boiling, nucleate boiling and single- 
phase convection that carries a close similarity to the pool boiling mechanisms. In general, the 
cooling rate of the tube wall is proportional to the mass flow rate of the quenching flow. In 
microgravity, there is no stratified flow and because the removal of the gravitational force will 
send liquid filaments to the central core and replace them by low thermal conductivity vapor that 
significantly reduces the heat transfer from the wall. Thus, the chilldown process is even longer 
in microgravity. 

6. Patents, Publications and Presentations 

6.1 Patents 

There was no patent application. 

6.2 Publications and Presentations 
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6. Development of Nanocrystalline Complex Metal Hydrides for 
Hydrogen Storage 

Task PI: Dr. Fereshteh Ebrahimi, Material Science & Engineering, University of Florida 
Graduate Student: Sankara Tatiparti, Material Science & Engineeringm, University of Florida 
Research Period: August 3, 2004 to March 31, 2007 

Abstract 

The objective of this work was to synthesize, via electrodeposition techniques, nanocrystalline 
Mg and Al alloys which can potentially be used for producing metal hydrides with low 
hydrogenation/dehydrogenation temperature, fast kinetics and reasonable gravimetric storage 
capacity. We have been successful in fabricating Al-Mg alloys in the form of porous powder. 
The composition of the deposits was varied with changing the amount of magnesium ions 
introduced into the electrolyte by a pre-electrodeposition process. The effects of current density 
and cathode material (copper versus graphite) were investigated. It was shown that the alloy 
powders, depending on the magnesium content, consisted of fee Al-Mg and/or hep Mg-AI 
supersaturated solid solutions. No intermetallics were found in the as-deposited powders. The 
maximum solubility of Mg in Al was found to be approximately 20at%, but supersaturated Mg 
phase with as high as 40at%AI was observed. TEM studies revealed that the grain size of these 
powders is ultra-fine. Low temperature heat treatments elucidated that the supersaturated 
aluminum has a higher thermal stability than the supersaturated magnesium alloy does. We 
have been able to secure funding from National Science Foundation (DMR-0605406) for 
continuation of this work. 


Accomplishments 


1. Electrodeposition set-up 

We finished the installation of the glove box and the electrodeposition cell was prepared to fit 
our needs. Figure 1 presents the electrodeposition set-up. Since Al and Mg cannot be deposited 
from aqueous solutions, we have to use organometallic-based electrolytes. These materials are 
highly sensitive to the oxygen and moisture, and therefore, an atmosphere of less than 5 ppm 
moisture and 1 ppm of oxygen has to be maintained during mixing of chemical and 
electrodeposition using the glove box. 

We employed a rotating cylinder cell design with pure copper as the substrate. The potential or 
current was controlled using a PAR 273 Potentiostat/Galvanostat, which was interfaced with a 
computer. For the reason of moisture and oxygen sensitivity of the electrolyte a non-aqueous 
Ag/AgCI reference electrode in a 1 molar LiCI in ethylene glycol monobutyl ether solution was 
used. The copper electrode was electropolished outside the glove box using a 82.5% 
phosphoric acid + 17.5% de-ionized water solution at a potential of 0.6V (SCE) at room 
temperature using a stainless steel cathode and then it was transferred inside. 
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2. Electrodeposition of pure aluminum 

Prior to the deposition of Al-Mg alloys, the electrodeposition of the aluminum was attempted to 
serve as trial run on the organometallic chemical to be used as the main component of the 
electrolyte. The electrolyte used was made of NaF, AIEt 3 , Toluene as suggested by Lehmkuhl 
et al. [1]. Electrolyte compositions with two different compositions, namely NaF+(2 or 4mol) 
AIEt 3 (where Et = -C 2 H 5 ) +(3 or 6mol) Toluene were studied. Toluene a solvent and NaF + AIEt 3 
form a complex according to the equation 1 . 

NaF + 2AI(Et) 3 Na[Et 3 AI-F-AIEt 3 ] (1 ) 



Figure 1 . Pictures of deposition set-up (a) and the glove box (b) used for the deposition of 
pure Al. 

Prior to electrodeposition, cyclic voltammetry was conducted in order to select the suitable 
potential for the deposition of aluminum powders. The cyclic voltammetry curve is shown in 
Figure 2. The fluctuations in the curve are due to the rotation of the copper cathode at 400 rpm. 
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Figure 2. Cyclic voltametery results with a scan rate of 3.33 V/s at 75°C and rotation speed of 
400 rpm. Potential is measured versus Ag/AgCI reference electrode. 

A potential of -1.5 V was selected for the electrodeposition of aluminum based on the 
voltametery results. The deposition was conducted potentiostatically at 75 °C for 15 minutes 
with cathode rotation of 400 rpm. Figure 3(a) presents the current density versus time curves for 
the two electrolyte compositions. Since only complexes liberate the Al +3 ions to be deposited, 
according to equation (1), the electrolyte containing the 1:2 ratio of Na:AI(Et) 3 will have a higher 
concentration of complexes than the electrolyte with 1:4 ratio. The higher current density 
observed using the former electrolyte can be attributed to the higher concentration of the 
available Al +3 ions. Figure 3(b) shows a low magnification picture of the copper electrode with 
the aluminum deposit on it. 



0 200 400 600 800 1000 1200 

Time (s) 



Figure 3. (a) Current density as a function of time for potentiostatic deposition of pure Al using 
two different electrolyte compositions, (b) Copper electrode with deposited aluminum. 
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The SEM (scanning electron microscope) images of the deposits of pure Al are shown in 
Figure 4. It can be observed that the morphology and the size of powders depend on the 
electrolyte composition. A higher concentration of available Al +3 (NaF: AI(Et) 3 ratio of 2) and the 
resultant higher current density led to a deposit with a much smaller powder size (Figure 4a). 
This observation is consistent with the theoretical prediction that the nucleation rate should 
increase with increasing the current density. 



Figure 4. SEM pictures of deposits made with (a) NaF: AI(Et) 3 ratio of 2, and (b) NaF: AI(Et) 3 
ratio of 4. 

Figure 5 presents the energy dispersive spectra (EDS) showing that the deposit consisted of 
pure aluminum. The x-ray diffraction (XRD) profile of the pure aluminum powder is shown in 
Figure 6. 
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Figure 5. EDS spectra of an Al deposit. The Cu peak is due to the residue from the substrate. 
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Figure 6. XRD profile for pure aluminum deposits. 


3. Electrodeposition of aluminum-magnesium alloys 

For making Al-Mg alloys, magnesium can be introduced as diethyl- or dimethyl-magnesium into 
the electrolyte [2], however, the transportation of these materials is not allowed in the US. An 
alternative method to introduce Mg is to use a magnesium anode, which dissolves during the 
deposition and provides the necessary Mg for forming complexes. It has been demonstrated 
that an electrolyte with Na[AIEt 4 ]/2Na[Et 3 AI-FI-AIEt3]/2AIEt3/6 toluene will allow the dissolution of 
Mg at anode [3,4], The Na[AIEt 4 ] can be produced by reaction of 3Na with 4AIEt 3 or by reaction 
of NaH with AIEt 3 in the presence of ethylene [4], In this case a pre-electrolysis should be 
carried out to increase the Mg content in the electrolyte before the actual deposition. 

We attempted to prepare several electrolytes using different procedures. After many trials we 
were successful in establishing the correct procedures for preparing the electrolyte and pre- 
electrolysis. We reduced the electrolytic cell size from 100 ml to 30 ml to decrease the chemical 
waste. Based on our studies, we improved our technique for producing Na[AIEt 4 ] and selected 
galvanostatic pre-electrolysis at a current density of 60 mA/cm 2 and temperature of 90 °C with a 
cathode rotation speed at 200 rpm. 
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Figure 7. (a) Low and (b) high magnification SEM pictures of the deposit with 5.5 at.%Mg. 
(c) EDS profile. 


In order to minimize the deposition of sodium we added an excess of 0.5 M of AIEt 3 to our 
electrolyte containing Na[AIEt 4 ], 2M Na[Et 3 AI-H-AIEt 3 ], 2M AIEt 3 and 6M toluene. We varied the 
duration of the pre-deposition from 30 min to 4 hours 25 min to establish the composition of the 
deposits as a function of the pre-electrolysis time. The composition as a function of the time is 
shown in Figure 7. It should be noted that this concentrations of Mg and Na were evaluated 
using EDS in SEM and therefore, they represent approximate values. The curve has a parabolic 
shape suggesting that a steady state can be reached at long times (or high Mg concentrations) 
where the amount of Mg dissolved from the anode becomes equal to the amount of Mg 
deposited, which results in a constant deposit composition. It is expected that this saturation 
concentration can be modified using anodes made of Al as well as Mg. The amount Na was 
reduced but it was not completely eliminated. 
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Figure 8. The approximate composition of the Al-Mg alloys as a function of the duration of 
electrodeposition. 



All three deposits whose compositions are given in Figure 8 exhibited dendritic growth as shown 
in Figure 9. The size of the dendrite arms became smaller and more directional as the amount 
of Mg in the alloy was increased. 

The XRD results for the Al-Mg deposits are presented in Figures 10 through 12. These results 
indicated that the 20 at.% Mg deposits was in a supersaturated condition containing only the 
solid solution of magnesium in aluminum. The 40 at.% Mg deposit consisted of supersaturated 
solid solution of magnesium in aluminum along with AI 3 Mg 2 intermetallic. The 60 at.% Mg 
deposit contained both AI 3 Mg 2 and AI 12 Mg 17 intermetallics. Again, the presence of Cu is due fact 
that some copper was scraped from the substrate during the removal of the deposited alloys. 

Our objective is to fabricate supersaturated nanocrystalline Al-Mg alloys with approximately of 
25-35 at.% Mg. A recent study of sputtered Al-Mg alloys indicate that this range of composition 
may result in the formation of magnesium alanate, Mg(AIH 4 ) 2 upon hydrogenation [5]. We have 
been successful in preparing supersaturated solid solution of Al-Mg. The next step is to design 
deposition schemes that result in more precise composition control with least variation as well 
as smaller powders with minimum grain size. 
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Figure 9. SEM pictures showing the morphology of the deposits with 
nominal compositions of (a) 20 at.%Mg, (b) 40 at.%Mg, and (c) 60 at.%Mg. 



Figure 10. The XRD profile for powders with an average composition 
of approximately 20 at.% Mg in Al. 
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Figure 1 1 . The XRD profile for powders with an average composition 
of approximately 40 at.% Mg in Al. 



Figure 12. The XRD profile for powders with an average composition 
of approximately 60 at.% Mg in Al. 

Milestones 

1 st - Research and design the electrodeposition set-up. 

2 nd - Order and prepare components. 

3 rd - Initial attempt for depositing Mg-AI alloys. Microstructural characterization of the deposits. 

4 th - Establish the relevant electrodeposition parameters for obtaining the desired 
nanocrystalline Mg-AI alloys. Characterize the properties of the alloys. 
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Products and Deliverables 


1 st - Glove box delivered and installed, interim report. 

2 nd Electrodeposition set-up completed, interim report. 

3 rd - Nanocrystalline Mg-AI material fabricated, interim report. 

4 th - - Properties of the Mg-AI alloy characterized and the electrodeposition parameters adjusted 
for improved properties, final report. 


References 

1 . H. Lehmkuhl, K. Mehler, U. Landau, “The principles and techniques of electrolytic aluminum 
deposition and dissolution in organoaluminum electrolytes”, H.Gerischer, C. Tobias, 
“Advances in electrochemical science and engineering”, 3, VCH, p. 163. 

2. Mayer A, “Electrodeposition of aluminum, aluminum/magnesium alloys, and magnesium 
from organometallic electrolytes”, Electro Chem. Soc. 137: 2806-2809 (1990). 

3. Lehmkuhl H, Mehler K, Reinhold B, Bongard H, Tesche B “Deposition of aluminum- 
magnesium alloys from electrolytes containing organo-aluminum complexes" Adv. Engg. 
Mater. 3:412-417(2001). 

4. Lehmkuhl H, Mehler K, Reinhold B, Bongard H, Tesche B “Deposition of aluminum- 
magnesium alloys from electrolytes containing organo-aluminum complexes” 
Materialwissenschaft und Werkstoffechnik 31 :889-898 (2000). 

5. R. Gremaud, A. Borgschulte, W. Lohstroh, H. Schreuders, A. Z' uttel, B. Dam, R. Griessen, 
Journal of Alloys and Compounds 2005, 404-406 (2005), pp. 775-778. 


NASA/CR— 2008-21 5440/PART3 


198 



7. Lithium Borohydride for Hydrogen Propellant Storage 

Task PI: Dr. Samim Anghaie, Nuclear & Radiological Engineering Department 
Graduate Student: Anne Charmeau, Nuclear & Radiological Engineering Department 
Undergraduate Student: Laura Padilla, Nuclear & Radiological Engineering Department 
The final report was submitted for this project. It is given again. 

Research Period: August 3, 2004 to April 1, 2006 
Projects Goals 

The proposed research task will be focused on the development and characterization of lithium 
borohydride for the dual hydrogen propellant storage and neutron shielding applications. Lithium 
borohydride uniquely combines two important properties that are pivotal to the operation of 
space propulsion systems; it has the highest hydrogen storage density and the highest neutron 
absorption cross section when fully enriched 6 Li and 10 B are used. 

Abstract 

Lithium borohydride features unique properties for hydrogen storage and radiation shielding. 
Room temperature properties of lithium borohydride are measured and compared with data 
available in open literature. Measurement of high temperature properties of lithium borohydride 
requires additional safety measures and operation under inert gas environment. Radiation 
shielding properties of lithium borohydride with and without the use of enriched lithium and 
boron isotopes were quantified and analyzed. Results indicated that significant improvement in 
neutron shielding effectiveness can be achieved by replacing natural lithium borohydride with 
fully enriched 6 Li 10 BH 4 . This is especially true at lower neutron energies where the neutron 
absorption cross-sections for 6 Li and 10 B are much higher than their naturally occurring isotopic 
forms. 


Introduction 

Storage of a gram of hydrogen gas at atmospheric pressure would require about 1 1 liter of 
space. So for convenience the hydrogen gas must be intensely pressurized to several hundred 
atmospheres or it must be stored in liquid form under cryogenic temperatures. These options 
are not practical for long duration space missions. A viable solution to these difficulties is 
storage of hydrogen in hydride form. Among all hydride materials, lithium borohydride (LiBH 4 ) 
has the highest hydrogen density and the lowest mass density in solid phase (0.82 g/cm 3 ). In 
addition, the unique nuclear properties of basic elements of the LiBH 4 would make it the most 
effective shielding materials for neutrons. Hydrogen is one of the best material for moderation 
and shielding of neutrons. Lithium and boron in their low isotopic forms, 6 Li and 10 B, are among 
the best neutron absorbing materials. The relative aboundance of 6 Li and 10 B are 7.5% and 
19.9%, respectively. The application of 6 Li 10 BH 4 could potentially bring about a revolutionary 
reduction in the size, weight, and the cost of future space propulsion systems for long duration 
manned missions to the solar system. 
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The focus and specific tasks of the research included: 

• Research Focus 

- Development and characterization of lithium borohydride for the dual hydrogen 
propellant storage and neutron shielding applications 

• Research Tasks 

- Measure and compile physical and thermodynamic property data for lithium 
borohydride at room temperature. 

- Perform thermal stability analysis of lithium borohydride at temperatures up to the 
total dissociation point (-700K) for the compound. 

- Performance and microstructural characterization of lithium borohydride under 
intense radiation field. 


Experimental Facility 

Due to high combustibility of lithium borohydride in presence of air and moisture, all testing and 
preliminary characterizations were conducted in the Innovative Nuclear Space Power & 
Propulsion (INSPI) see-through low vacuum chamber. X-ray Diffraction (XRD - Philips APD 
3720) technique was used for identification of chemical compound. This device uses powder 
diffraction standards files of d-spacings and intensities of known compounds to determine 
crystallographic information such as crystal structure, reflections and lattice parameters. 

Accomplishments 

An extensive literature search was conducted to update the existing database on lithium 
borohydride. The crystalline structure of LiBH 4 has been studied using synchrotron soft x-ray 
diffraction (J-P Soulie et.al. J. Alloys Compnd., 346, 200 (2002)). At ambient conditions LiBH 4 
has orthorhombic symmetry. At higher temperatures (T > 100C), LiBH 4 undergoes phase 
transition to hexagonal structure. In both structures, four hydrogen atoms are arranged around 
boron in a tetrahedral configuration. 

First principal calculations for prediction of fundamental properties of LiBH4 have been 
performed by K. Miwa et. al. (Phys. Rev. B, 69, 245120 (2004)) using the ultrasoft 
psedopotential method based on density function theory. The generalized gradient 
approximation approach has been used for the exchange-correlation energy and prediction of 
LiBH 4 fundamental properties such as the structural properties, electronic properties, dielectric 
properties, vibrational properties, and heat of formation. The calculation and analysis suggested 
a nearly ideal tetrahedral configuration of BH 4 , which is in contrast to recently measured 
structural properties that suggest strongly distorted configuration. Calculation of electronic 
structure indicated that Li atoms are ionized as Li+ cations. 

Known properties of lithium Borohydride pertinent to hydrogen storage include: 

- 19.6% H content 

- Molecular weight 21 .78 g/mol (LiBH 4 ) 

20.22 g/mol (6Li10BH 4 ) 

- Density 0.66 g/cm 3 

- Melting point @ 0.1 MPa 541 K 

- Boiling point @ 0.1 MPa 551 K 

- Decomposition point @0.1 MPa 473 K 

- White crystalline solid 
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- X-ray power diffraction: Orthorombic symmetry 

- Stable in dry air or vacuum (T < 473 K) 

- Transition to hexagonal structure ~381K 

- Hydrogen desorption peak 700K 

- Primary solvent Tetrahydrofuran (THF) 

The primary methods of preparing high purity LiBH4 is through the following chemical reactions: 
NaBH 4 + LiB LiBH 4 + NaB ► 

(Under nitrogen environment, 18 h at 298 K and .5 h at 308 K) 

Li + B + 2H 2 LiBH 4 ► 

(Heat of formation ~ - 194 KJ/mol) 

The hydrogen desorption reaction proceeds at temperatures above 573 K with Si0 2 as catalyst: 
LiBH 4 LiH + B + 3/2 H 2 ► 

For the purpose of this project, the destabilization of LiBH 4 to decrease the dehydriding 
temperature is a primary focus of research. Based on the published data, it is suggested that 
the charge compensation by Li+ cations play a pivotal role in the stability of the internal bonding 
of [BH4]- anions. It is also expected that the suppression of the charge transfer by the partial 
substitution of Li by other elements will be effective to lower the dehydriding temperature. 

Thermal stability experiment was conducted to measure the melting point and stability of natural 
LiBH 4 at temperatures ranging from 300 K to 800 K. 500 g of ultrahigh pure natural LiBH 4 in 
powder form was acquired. Due to high flammability of LiBH 4 , small packs of 2 g powder were 
prepared to minimize potential damage to the inert environment furnace system. During the 
preparation of test samples, it was observed that even at very low humidity the rate of moisture 
absorption is too high. LiBH 4 is soluble in water at room temperature. During weighting of LiBH 4 
under normal atmospheric conditions, the air moisture dissolved LiBH 4 particles to form small 
droplets. The process was moved to an inert gas chamber to minimize reaction with moisture. 
X-ray analysis was performed to obtain diffraction pattern for both pure and dissolved LiBH 4 
samples. Figure 1 shows the x-ray diffraction pattern for LiBH 4 samples prepared under inert 
gas atmosphere and at room temperature (~20 °C). 
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Figure 1. LiBH4 x-ray diffraction pattern. 

High temperature measurement of the LiBH 4 properties was hindered by the safety concern 
associated with high combustibility in presence of air and moisture. In two occasions during 
handling of small LiBH 4 samples from the inert gas chamber to the analytical lab and vacuum 
furnace, spontaneous combustion resulted in small and limited lab fire. Though none of the lab 
fires resulted in any property damage or personal injury, the high temperature experiment with 
LiBH 4 was suspended pending the review of safety procedures. Further analysis of the 
experimental procedure and properties of LiBH 4 yielded the conclusion that high temperature 
testing only should be performed under inert gas environment. This change would have 
necessitated moving the furnace and basic analytical instruments under the controlled 
environment of an inert gas chamber. Due to significantly additional costs, continuation of this 
task was postponed to the nest funding cycle of this project. 

In addition to hydrogen storage, LiBH 4 could be used for shielding of space nuclear reactors. To 
assess and evaluate the shielding capability of LiBH 4 , a computational study was conducted 
using the industry standard MCNP-5 computer code system [3] . In particular, cross-sections for 
natural lithium and boron in LiBH 4 were replaced with highly enriched 6 Li and 10 B, respectively. 
The neutron absorption cross-sections for these isotopes are much higher than the other 
naturally occurring isotopes of these elements. This is a key factor in reducing the overall weight 
of shielding materials that are needed for protection against radiation originated from external or 
internal sources. 

A cylindrical model with a neutron emitting line at the center was developed to evaluate and 
compare shielding properties of LiBH 4 with and without isotopic enrichment. Two different set of 
runs were conducted using this geometry, one for natural LiBH 4 and one for 6 Li 10 BH 4 . Both 
materials were modeled under the same conditions. In order to check the effect of the neutron 
energy and depth of the material on its performance as a shield, values of both the cylindrical 
model radius and the neutron source energy were varied. The radius of the cylinder varied from 
10 to 50 cm that covers a broad range of proposed space nuclear reactor designs. The energy 
of the linear neutron source was set at four different points, 0.010, 0.100, 0.500 and 1 MeV, 
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which in the lower end of the fission spectrum that is shown in Figure 2. The height of the 
cylinder was maintained constant at 20 cm. Results of the analysis indicated that though both 
natural and enriched lithium borohydride materials have excellent shielding capabilities, 6 Li 10 BH 4 , 
has significantly higher shielding capability for low energy neutrons. 



Figure 2. Maxwellian distribution of fission neutron energy (fission spectrum) 

Results of the analysis indicated that though both natural and enriched lithium borohydride 
materials have excellent shielding capabilities, Li 6 B 10 H 4 , has significantly higher shielding 
capability for low energy neutrons. This is expected because Li and B neutron absorption cross- 
sections increase with decreasing neutron energy (especially for Li 6 and B 10 ). However, the 
shield still proves to be highly effective in higher energies. At 1 MeV, according to our 
calculations, only 0.152 % of the neutrons make it through a 10 cm shield of natural LiBH 4 . This 
is impressive since 10 cm is the thinnest shield we modeled (therefore the most ineffective one) 
and natural LiBH 4 is not as good a shield as Li 6 B 10 H 4 This value drops to 0.109% when using 
Li 6 B 10 H 4 . 

Figure 3 shows the calculated value of neutron flux penetrated through 40 cm of lithium 
borohydride with and without isotopic enrichment at different energies. The energy integrated 
values of the neutron flux penetrated through 10 cm of lithium borohydride is shown in Figure 4. 
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Figure 3. Neutron flux and energy penetrated through 40 cm of lithium 
borohydride with and without isotopic enrichment. 
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Figure 4. Neutron flux penetrating through 10 cm of lithium borohydride as a 
function of source neutron energy. 
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Further analysis of results indicated that as the radius of the cylinder increases the relative 
effectiveness of the shield will also increase. The amount of neutrons transmitted through a 
20 cm shield is 25% of what is transmitted through a 10 cm shield. The amount of neutrons 
transmitted through a 50 cm shield is 16% of that transmitted through a 20 cm shield. These 
percentages remain fairly constant with increasing energy and they practically remain 
unchanged going from natural to enriched LiBH 4 . Another trend shown in these results is the 
significant improvement in shielding effectiveness by replacing natural lithium borohydride with 
fully enriched 6 Li 10 BH 4 . This is especially true at lower energies where the neutron absorption 
cross-sections for fully enriched 6 Li and 10 B are much higher than their naturally occurring 
isotopic forms. For 100 keV source neutrons, 50% fewer neutrons penetrate through the shield 
when natural LiBH 4 is replaced with fully enriched 6 Li 10 BH 4 . For 1 MeV source neutrons, the 
percentage reduction in penetrated flux ix approximately 30%. 

In summary, results of this study on hydrogen storage and shielding capabilities of lithium 
hydroboride have yielded the following conclusions; 

1 . Lithium Borohydride for Hydrogen Propellant Storage: 

• Lithium Borohydride features superior properties for hydrogen storage, in particular, for 
once use and non-recycled application of hydrogen. 

• Chemical stability and high energy reaction with water and air limits high temperature 
and terrestrial applications of Lithium Borohydride. 

• Measured properties at room temperature agreed with reported data. 

• X-ray analysis was performed to measure the moisture absorption and reaction. 

• White, odorless powder, water-reactive and is corrosive in the presence of moisture. 

• Reaction with water generates flammable hydrogen gas, causing ignition or explosion. 

• Thermal decomposition of Lithium Borohydride produces hydrogen gas, lithium and 
boron oxides, and borane. 

• Lithium Borohydride will react with water, carbon dioxide, and oxygen in air to form other 
inorganic compounds. 

• High combustibility of Lithium Borohydride in presence of air and moisture necessitates 
additional safety measures and procurement of new equipments that have not been 
foreseen at the outset of this project. 

2. Unique radiation Shielding Features of Lithium Borohydride 

• LiBH 4 with natural isotopic constituents is an excellent neutron shielding material. 

• Enrichment of natural B and Li to 10 B and 6 Li, respectively, further enhances neutron 
shielding property of LiBH 4 

• Monte Carlo analysis is performed to quantify the enhanced neutron shielding behavior 
of 6 Li 10 BH 4 . 
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Work in progress includes effect of radiation - neutron and gamma-rays - on properties 
of Lithium Borohydride 


3. Potential Applications of Lithium Borohydride: 

• Compact hydrogen storage. 

• Light weight & compact radiation shielding for space reactors. 

• Dual use as propellant storage for Nuclear Thermal Rocket (NTR) and rocket fuel for lift 
up from Moon or Mars surface. 
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